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I.  INTRODUCTION 

The  concept  of  observability  introduced  by  Kalman  in  1960  [Ref.  1],  is,  along 
with  the  dual  concept  of  controllability,  one  of  the  most  significant  ideas  to  come  out 
of  the  state  space  approach  to  control  theory.  The  concept  of  observability  presented 
in  most  of  the  literature  available  today,  gives  only  an  answer  to  whether  the  system  is 
observable  or  not.  However,  another  question  was  raised  by  Brown  in  1966  [Ref.  2], 
opening  the  issue  of  how  observable  a  given  system  is.  Are  all  the  states  equally 
observable?  If  not,  what  are  the  relative  degrees  of  observability  of  the  state  variables? 
Since  1966,  this  issue  has  been  investigated  in  [Ref.  3],  [Ref.  4],  [Ref.  5]  and  [Ref.  6] 
among  others.  In  the  majority  of  theses  investigations  the  goal  has  been  either  the 
state-space  reduction  of  the  model  or  frequency  domain  model  reduction.  Silverman 
has  probably  been  one  of  the  most  active  researchers  on  this  subject,  using,  in  several 
of  his  papers,  the  technique  of  rotating  the  coordinate  system  in  such  a  way  that  a  so 
called  "balanced  realization"  is  obtained.  By  studying  this  new  realization  the  "nearly 
redundant"  states  can  be  distinguished  from  the  most  important  ones,  providing  the 
required  information  to  reduce  the  order  of  the  model. 

In  this  thesis  the  problem  of  how  observable  the  states  are  is  again  studied  using 
techniques  different  from  the  ones  mentioned  above.  Two  graphical  methods  are 
developed:  the  first  one  is  relatively  simple,  but  rather  helpful  in  understanding  the 
concepts  involved  in  ranking  the  state  variables  according  to  their  relative  degrees  of 
observability.  The  second  graphical  method,  based  on  some  ideas  given  in  [Ref.  7],  is 
more  elaborate.  Both  graphical  methods  are  compared  against  a  third  method 
described  in  [Ref.  7]  whose  development  is  explained  in  this  thesis. 

The  tests  to  find  the  relative  observability  for  each  state  variable  provide  the 
information  to  develop  a  method  that  enables  the  improvement  of  the  degree  of 
observability  of  the  least  observable  state.  Although  all  the  basic  derivations  assume  a 
single  output  regulator  problem,  the  generalization  to  multi-output  non-regulator 
systems  is  also  made.  This  thesis  is.  therefore,  divided  as  follows:  in  Chapter  II  a  brief 
and  general  definition  of  observability  is  given  followed  by  two  single  output  numerical 
examples.  Chapter  III  begins  with  a  definition  of  relative  observability  followed  by  the 
presentation  of  methods  to  evaluate  the  relative  degree  of  observability  of  a  system. 


Three  algorithms  are  presented  leading  to  a  Fortran  computer  program  whose  listing  is 
in  Appendix  A.  In  Chapter  IV  an  algorithm  to  improve  the  relative  degree  of 
observability  of  a  weakly  observable  state  is  developed  for  the  single  output  regulator 
problem.  This  algorithm  leads  to  a  Fortran  program  that  is  contained  in  Appendix  B. 
The  generalization  for  multi-output  servo  systems  is  also  discussed.  Finally  in  Chapter 
V  the  performances  of  the  three  methods  of  Chapter  III  are  compared  as  well  as  the 
performance  of  the  algorithm  developed  in  Chapter  IV.  Some  particular  examples  used 
with  the  computer  programs  are  presented  and  analized  in  some  detail. 


II.  THE  CONCEPT  OF  OBSERVABILITY 

A.  DEFINITION 

A  system  is  said  to  be  completely  observable  if,  from  the  measurements  of  the 
output  y(t)  over  a  finite  period  of  time,  t  e[0,t()].  every  initial  state  can  be  determined 
[Ref.  8:  p.  764],  [Ref.  9:  p.  107],  In  other  words,  if  every  transition  of  the  states  affects 
every  element  of  the  output  vector,  the  system  is  completely  observable  [Ref.  8:  p.  764], 

B.  THE  STATE-OBSERVABILITY  PROBLEM 

Consider  a  linear  time-invariant  system  described  by  equations  2.1  and  2.2 

x(t)  =  Ax(t)  +  Bu(t)  (2.1) 


y(t)=Cx(t)  (2.2) 

where: 

\(t)  is  a  nxl  column  vector  of  the  state-variables. 

x(t)  is  a  nxl  column  vector  of  the  time  derivatives  of  the  state  variables 

A  is  a  nxn  matrix. 

B  is  a  nxp  matrix 

u  is  a  pxl  column  vector  of  the  inputs  to  the  system. 

y(t)  is  a  mxl  column  vector  of  the  output  variables  of  the  system 

C  is  a  mxn  matrix 

The  effect  of  the  input  signal  u(t)  can  always  be  determined,  if  x(0)  can  be 
computed  from  the  measurements  of  the  output,  i.e.,  if  the  system  is  completely 
observable.  Therefore,  there  is  no  loss  of  generality  by  assuming  that  u(t)  =  0  [Ref.  9:  p. 
107].  Under  this  assumption,  to  evaluate  a  system's  observability  a  matrix  Q.  called 
the  observability  matrix,  must  be  formed.   [Ref.  10:  p.  108]. 

Q=[CT    ATCT    (AT):CT     ...   (A1)11"^1]  (2.3) 
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Each  of  the  submatrices  separated  by  spaces  in  equation  2.3  is  nxl.  Therefore. 
there  will  be  n  columns  in  the  Q  matrix;  this  matrix  will  be  square  (nxn)  when  there  is 
only  a  single  output.  If  the  n  columns  of  Q  are  all  linearly  independent,  i.e..  if  the  nxn 
matrix  Q  is  nonsingular  the  initial  state  vector  x(0)  can  always  be  found  from  the 
measurements  of  the  output  [Ref.  10:  p.  105].  On  the  other  hand,  if  Q  is  singular,  (not 
of  full  rank),  some  columns  of  Q  will  be  linearly  dependent  on  others,  and  only  special 
choices  of  y(t)  will  produce  the  required  solution  [Ref.  10:  p.  105]. 

C.       EXAMPLE 

1.  Observable  system 

Consider  the  single  output  linear  time  invariant  system  described  by: 


Where: 


x(t)  = 

Ax(t)+Bu( 

t) 

y=Cx(t) 

1     2 

0 

0 

c 

A  = 

2     3 

0 

B  = 

0 

0     1 

1 

1 

Q  is  equal  to 

- 

1 

1     7 

Q  - 

0 

3    12 

1 

1     1 

-[io    .] 


This   matrix   Q   is   nonsingular,    i.e..    of  full   rank,   therefore   the    system  is 
observable. 

2.  Non-observable  system 

Consider,  again  the  same  state  representation  where  the  matrices  A,  B  and  C 
are  as  follows: 


A  = 


0     2 

0 

0    0 

-2 

0     2 

1 

B  = 


C  =  [o    0    1 ] 
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The  matrix  Q  is  equal  to: 


Q 


0    0    0 

0  2     2 

1  1    -3 


In  this  case  the  matrix  Q  is  singular,  i.e.,  its  determinant  is  equal  to  zero.  Its 
rank  is  also  less  than  three,  (actually  two),  because  the  maximum  submatrix  of  Q  with 
determinant  different  from  zero  is  a  2x2  matrix.  Therefore,  in  this  case  there  is  one 
state  that  is  not  observable.  The  identification  of  that  non  observable  state  is  not 
possible  just  by  performing  the  previously  described  test,  where  the  answer  "yes-no"  is 
given,  with  no  indication  of  how  observable  the  system  is  and  which  are  the  most 
observable  states.  In  the  following  Chapters  methods  will  be  discussed  and  presented 
where  answers  to  these  questions  will  be  given. 


III.  THE  CONCEPT  OF  RELATIVE  OBSERVABILITY. 

A.  DEFINITION 

The  problem  of  relative  observability  deals  with  finding  an  algorithm  or  figure  of 
merit  for  each  state-variable,  which  will  reflect  how  observable  the  state-variable  is. 

Measuring  y(t)  at  times  t  =  ti,  t-> tn,  one  wants  to  be  able  to  evaluate  which  state 

components  can  be  most  accurately  determined  from  the  n  set  of  measurements  of  the 
output.  In  other  words,  how  an  error  in  the  measurements  is  reflected  in  determining 
the  initial  states?  If  a  small  error  in  the  measurements  produces  a  large  error  in  the 
calculated  state,  this  state  is.  therefore,  difficult  to  determine  from  the  measurement  o[ 
the  output.  The  output  has  a  "fuzzy  image"  of  the  state  that  is  "hidden"  by  the  large 
error.  If  a  small  error  in  the  measurements  produces  a  small  error  in  the  evaluation  of 
a  particular  state,  this  state  will  be  evaluated  more  accurately  than  the  others  from  the 
measurements  of  the  output.  Due  to  the  small  error  the  "image"  of  the  state  at  the 
output  is  less  "fuzzy". 

B.  METHODS  FOR  EVALUATING  RELATIVE  OBSERVABILITY 

The  three  methods  that  will  be  described  for  evaluating  the  relative  degree  of 
observability  were  developed  based  on  the  concept  given  previously:  the  projection  of  a 
measurement  error  in  computing  the  initial  state. 

Starting  with  equations  2.1  and  2.2  and  assuming  the  input  to  be  zero,  equation 
2.2  is  differentiated  n-1  times  and  equation  2.1  is  substituted  as  shown  below 
[Ref.  10:  p.  106]. 

y(t0)=Cx(t0)  (3.1) 

y(t0)=Cx(t0)=CAx(t0) 

y(t0)=CAx(t0)=CA2x(t0) 


CT/t0)=CAn-1x(t0) 
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The  result  is  n  equations  with  n  unknowns.   It  can  be  seen  that  the  coefficient 

vectors   C,   CA,   CA  , CA'       '   are   row  vectors   with   the   same   elements   as   the 

respective  column  vectors  C    ,  A'C    .(.VrC1 (A1)11"^1   of  the  Q  matrix.  Thus, 

the  elements  of  the  columns  of  the  Q  matrix  are  the  coefficients  of  the  equations 
relating  the  initial  states  with  the  measured  quantities.  Equation  3.1  can  be  written  as 

yd(t())  =  QTx(t0)  (3.2) 

where  y^n)  is  the  column  vector  of  the  output  and  its  n-1  derivatives  evaluated  at  t  = 

The  fact  that  the  output  equation  2.2  must  be  differentiated  n-1  times  to  obtain  n 
equations  and  n  unknowns,  implies  that  the  output  has  to  be  measured  n  times.  If  the 

system  has  only  one  output,  as  in  this  case,  the  state  vector,  x(tQ),  can  be  found  by 

T  T 

inverting  the  Q     matrix.  The  inverse  of  Q     only  exists  if  its  determinant  is  nonzero  or, 

in  other  words,  if  the  rank  ol"  the  matrix  is  equal  to  its  order.  However,  if  the 
coefficients  of  any  pair  of  equations  are  approximately  proportional,  the  determinant 

T 

of  Q  will  be  a  small  number  and,  from  a  practical  point  of  view,  it  will  be  difficult  to 
invert  the  matrix  and  solve  the  system  of  equations.  A  matrix  that  has  a  small  number 
for  its  determinant  is  an  "ill-conditioned"  matrix.  Since  a  zero  value  for  the  determinant 
occurs  only  for  a  singular  matrix,  it  can  be  said  that  an  ill-conditioned  matrix  is  nearly 
singular. 

The  solutions  of  ill  conditioned  systems  are  very  sensitive  to  small  changes  in 
either  the  coefficients  of  the  matrix  or  the  values  of  the  forcing  functions.  Therefore,  if 

T 

the  rows  of  Q  are  nearly  proportional,  a  small  error  in  in  the  measurement  of  y^tQ), 
will  produce  a  large  error  in  x(tQ)  due  to  the  high  sensitivity  of  ill-conditioned  systems. 

C.       UPPER-BOUND  ERROR 

1.  Theoretical  background 

T 

The  degree  of  independency  of  the  rows  of  Q     can  be  used  to  quantify  the 

degree  of  observability  for  each  state  variable  [Ref.  7:  p.  12].    If  two  or  more  of  the 

T  T     1 

rows  of  Q     are  almost  linearly  dependent  the  entries  of  (Q    )     will  be  large  numbers. 

In  defining  the  degree  of  observability  for  each  state  variable  using  the  upper-bound 
error,  the  reciprocal  of  the  increase  in  the  error  of  the  calculated  state-variable  over  the 
observation  error  will  be  used  [Ref.  7:  p.   17].    The  reciprocal  will  be  used  because  a 
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small    number    associated    with    a    low   degree    of  observability    is    more    intuitively 
appealing  than  a  large  number. 

T 

Referring  to  matrix  equation  3.2,  normalization  of  the  rows  of  Q     will  be 
performed.    The  goal  of  this  operation  is  to  obtain  a  maximum  value  for  the  degree  of 

observability  equal  to  one.    To  accomplish  the  normalization  the  elements  of  yj<t())  are 

T 

divided   by   the   length   of   the   corresponding   row  vector   of  the   Q      matrix.     The 

T 
normalized  vector  is  defined  as  yn(tQ)  and  the  normalized  Q     as  P.    Equation  3.2  can 

now  be  written  as  : 


yn(t0)=Px(t0)  (3.3) 

The  vector  yn(tn)  consists  of  the  actual  value  of  the  measurements  ya(t(-j)  and 
an  error  e    .  thus,  equation  3.3  can  be  written  as  : 

ya(to)+em=Px(to)  {}A) 

Solving  equation  3.4  for  x(t/\)  by  taking  the  inverse  of  P  results  in 

x(tQ)=P-1va(t0»-P-1em  (3.5) 

Where: 

x(tn)  is  the  nxl  vector  of  the  computed  states  at  time  t(-,. 

P    ya(tQ)  is  the  nxl  vector  of  the  true  states. 

P"*e     is  the  nxl  vector  of  the  error  in  the  computed  states. 

Defining 

^'0»=P"1>a<V  (3'6> 

and  substituting  this  in  equation  3.5  yields 

K(t0)  =  xa(t0)+P-Iem  (3.7) 


or 
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x(t0)-xa(t0)=P"1em  (3.8) 

The  left  hand  side  of  equation  3.8  is  actually  the  error  in  the  computed  state, 
ec.  Substituting  this  notation   in  equation  3.8  gives  : 

ec=P-'em  (3.9) 

Recalling  what  has  been  said  about  the  dependency  among  the  rows  of  P,  it 
was  stated  that  if  the  dependency  was  high  the  system  was  ill-conditioned  and  the 
outcome  was  large  numbers  for  the  rows  of  P  .  Equation  3.9  shows  that  the  error  in 
the  computed  states  is  a  linear  combination  of  the  measurement  errors.  For  equal 
values  of  error  in  the  measurements,  i.e.,  ejm=  e->m  =  ...  =  enm,  the  entries  in  the  e 
vector  will  be  different  according  to  the  magnitudes  of  the  entries  in  the  corresponding 
rows  of  P.  The  larger  the  sum  of  magnitudes  of  the  entries  in  a  row,  the  larger  will  an 
error  in  y^n)  ^e  reflected  in  x(tQ). 

To  recapitulate  .  the  upper-bound  error  criterion  was  defined  as  the  reciprocal 
of  the  ratio  of  the  error  in  the  calculated   state-variable  to  the  observation  error. 


e 
'obsi 


§nh<i  = =  —  (3-10) 


e-  e- 

1C  1C 

eim 


Where: 

^obsi  *s  t^ie  degree  of  observability  of  the  ith  state  variable. 

e'm  is  the  error  in  the  ith  measurement. 

eic  is  the  error  in  computing  the  ith  state. 


Equation  3.9  can  be  rewritten  as  equation  3.11,  where  the  a-  are  the  entries  of  P    . 

elc  =  allelm  +  C£12e2m-f-4-alnenm  (3J1) 

e2c  =  a21elm+a22e2m+-"fa2nenm 


enc~  anlelm+  an2e2m+  -  +  annenm 
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Substituting  equation  3.11  in  3.10  yields 


6^-  = f>m 


,    •  = an        (3  1  "M 

(«ilelm  +  «i2e2m  + +  «inenm> 

The  main  goal  is  to  obtain  a  measure  of  the  magnitude  of  the  error  in  the 
computed  states,  for  an  equal  error  made  while  measuring  the  output.  Therefore,  it  is 
reasonable  to  consider  as  equal  all  of  the  e  .  If  one  does  not  do  so,  misleading  values 
for  dobs  for  the  state-variables  may  be  obtained,  because  the  comparison  between 

^obsT  ^obs2 ^obsn*  w^  ^e  ^one  unc^er  different  conditions. 

Consider,  then,  that  all  the  errors  in  the  measurements  are  equal.  Factoring 
out  the  em  in  the  denominator  of  equation  3.12  and  cancelling  this  value  with  the 
numerator,  yields  equation  3.13  : 

Sobsr-rT^-)  <3-13) 

— v  ij; 

Since  one  is  looking  for  the  maximum  possible  value  for  the  error  in 
computing  the  ith  state,  the  signs  of  the  elements  in  the  em  vector  must  be  such  that 
an  upper-bound  for  the  computed  error  e-  will  be  obtained.  Since  the  computed  error 
is  a  linear  combination  of  the  observation  errors,  the  upper-bound  for  that  error  will  be 
given  by  the  sum  of  the  absolute  values  of  the  row  entries  of  P    .  This  is  expressed  by 


8obsi=TTrT  (114) 

ZJuijl 

which  is  the  final  equation  for  evaluating  the  degree  of  observability  for  each  state- 
variable  using  the  upper-bound  error  method. 
2.  Example 

Consider  the  single  output  linear  time-invariant  system  described  by: 

x(t)  =  Ax(t)+Bu(t) 

y(t)  =  Cx(t) 
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Where: 

A  = 

"  o  r 

B  = 

r0 

_-4     -5 
STEP  ONE  -  Form  the  Q 

matrix: 

1 

c  = 


Q=       [CTATCT]      => 


J 


Q  = 


[0      1] 


0  -4 

1  -5 


STEP  TWO  -  Form  Q1  and  normalize  its  rows,  producing  the  matrix  P: 


QT=r  o   r 

_-4      -5 

STEP  THREE-  Find  P'1: 


=  > 


P  = 


0         1 
-0.625  -0.781 


PM=    -1.249-1.6 
1        0 
STEP  FOUR  -  Find  the  degree  of  observability  using  equation  3.13: 


1 


'obsl 


'obs2 


|-1. 249|  +  |-1. 6| 
1 


=  0.351 


=   1 


1       4-     0 

The  upper-bound  error  criterion  for  Xi  says  that  the  error  on  computing  xWtn)  will  be 
2.S5  times  bigger  than  the  errors  made  while  measuring  v(tn)  and  its  derivative  if  the 
assumption  that  both  measurement  errors  are  equal  holds.  The  upper-bound  error 
method  criterion  for  X2  tells  the  user  that  the  error  on  computing  x->(tn)  will  be  the 
same  as  the  one  made  on  the  measurements.  This  result  is  very  satisfying  because, 
looking  at  the  C  matrix,  one  can  see  that  the  state  x->  is  measured  directly  while  Xj 
influences  the  output  measurement  indirectly. 

D.       STANDARD  DEVIATION  ERROR  CRITERION 

1.  Theoretical  background 

Assume  now  that  the  errors  in  the  output  measurements  that  are  being  made 
are  statistically  independent  random  variables  whose  second  order  statistics  are  known. 
It  is  known  that  the  variance  of  a  linear  combination  of  statistically  independent 
random  variables  is  the  sum  of  the  coefficients  squared,  multiplied  by  the  respective 
variances  of  each  random  variable  [Ref.  11:  p.  100].    Thus,  taking  the  variance  of  the 
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ith  measurement  and  dividing  it  by  the  variance  of  the  error  in  computing  the  ith  state- 
variable  gives 


8obsi" 


Var[e: 


mJ 


Vartailelm+ai2e2m  + +  ainennJ 


(3.15) 


If  the  stated  assumption  holds,  i.e.,  if  the  measurements  of  the  output  are 
statistically  independent  random  variables,  the  following  equation  is  obtained: 


Var[e5m] 


'obsi' 


imJ 


ail"  Var[elml  +  ai2'  Var£e2nJ  + +  ain"  Var[ennJ 


(3.16) 


Here  again  the  assumption  of  equal  variance  is  made  because  this  is  the  only  way  of 
comparing  the  n  degrees  of  observability  under  the  same  conditions.  Factoring  out  the 
variances  in  equation  3.16  and  cancelling  with  the  numerator  results  in  equation  3.17 


1 
>obsi "     au2  +al22  + +  ain2 


(3.17) 


1 


obsi      £|0  ,2 


(3.18; 


Equation  3.18  is  the  final  equation  for  evaluating  the  degree  of  observability 
for  each  state-variable  using  the  variance  as  a  measure  of  the  error. 
2.  Example 

Starting  in  step  three  of  the  previous  example: 


p-U 


■1.249  -1.6 
1        0 


Find  the  degree  of  observability  using  equation  3.18: 

1 


8obsl~ 


*obs2 


1.249)-  +  (-1.6) 

1 


(D-  +    0 


2=0.243 
=    1 
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The  variance  criterion  for  evaluating  the  degree  of  observability  shows  that 
the  error  in  computing  Xj(t(j)  has  a  variance  4.12  times  bigger  than  the  one 
encountered  in  the  measurements.  The  variance  of  the  error  in  evaluating  x->(t())  is 
equal  to  the  variance  of  the  error  in  the  measurements. 

E.        GRAPHICAL  METHOD 

Another  interpretation  of  ill  conditioning  is  obtained  by  examining  graphically  a 
system  of  two  equations  and  two  unknowns.    Consider  the  following  system: 

QTx(t0)  =  y(t0)  (3.19) 


where: 

*"  1.01    0.99 


QT 


0.99    1.01 


T 

Substituting  the  matrix  Q     in  equation  3.19  yields 


1.01x1(t0)  +  0.99x2(t0)  =  y1(t0)  (3.20) 


0.99x1(t0)+  l.01x2(t0)  =  y2(t0)  (3.21) 

Equations  3.20  and  3.21  can  be  seen  as  the  equations  of  two  lines  with  the 
abcissa  Xj(tQ).  the  ordinate  x->(tQ)  and  the  x2(tQ)  intercepts  at  Xi(tQ)  =  0  being 
proportional  to  y'j(tQ)  and  y2(tQ)  respectively.  The  measurements  of  the  output  are 
Vj(tQ)  and  y2(tQ)  and  the  main  interest  is  to  evaluate  how  much  an  error  or  disturbance 
made  in  the  measurements  (em)  is  "amplified"  when  the  computation  of  the  initial 
states  takes  place.  In  Figure  3.1  the  lines  corresponding  to  equations  3.20  and  3.21  are 
presented.  The  lines  intersect  in  the  point  Xia(t/\),  x->Jtn)  which  is  the  actual  solution 
of  the  system  of  equations.  Suppose  that  some  errors  are  made  while  measuring  the 
outputs  yj(tQ)  and  y2(t^).  How  is  this  fact  going  to  change  the  position  of  the  original 
lines?  How  is  the  solution  of  the  original  system  going  to  be  perturbed?  As  in  the 
development  and  discussion  of  Ablin's  method  [Ref.  7],  the  errors  in  the  measurements 
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will  be  considered  as  equal  so  that  the  comparison  between  the  different  degrees  of 
observability  will  be  done  under  the  same  conditions.  Consider,  then,  that  the  value 
measured  for  yj(tQ)  exceeds  its  real  value  by  an  amount  that  is  equal  to  em<  In  this 
case  the  line  corresponding  to  equation  3.20  is  displaced  in  the  direction  of  positive 
x->Uq)  and  Xj(tQ)  to  a  position  parallel  to  its  original  one.  Conversely,  if  the  measured 
value  is  less  than  its  real  one  by  an  amount  that  is  equal  to  e_,  the  line  corresponding 
to  equation  3.20  is  displaced  to  a  position  parallel  to  its  original  one  but  now  in  the 
direction  of  the  negative  x^^q)  and  Xj(tQ).  These  two  lines  form  an  uncertainty  region 
around  the  line  corresponding  to  equation  3.20.  The  same  thing  occurs  for  the  line 
corresponding  to  equation  3.21  and  a  second  uncertainty  region  is  obtained.  As  can  be 
seen  in  Figure  3.1  the  intersection  of  these  two  regions  forms  a  parallelogram  shaped 
region  that  contains  all  the  possible  solutions  of  the  system  of  equations  if  the 
measurements  can  van-  between  >" itto)-1-  em  am*  ^  l^O^m  *or  ^e  *~irst  measurement  °f 
the  output  and  y^o)"1"6!!!  an^  ^^O^m  ^or  ^ie  secon(^  one-  The  boundaries  of  the 
parallelogram  shaped  region  contain  the  solution  of  the  system  of  equations  for  the 
particular  cases  where  the  errors  in  the  measurements  are  maximum  in  magnitude.  In 
Figure  3.1  eight  different  points  are  marked  on  the  boundary  region.  If  the  errors  in  the 
measurements  are  such  that  yi(t«)  is  measured  at  its  minimum  allowable  bound. 
namely  yi(tQ)-em  and  y->(tf))  is  measured  at  its  maximum  allowable  bound,  point  1  is 
obtained.  In  the  case  of  point  2  no  error  was  assumed  in  the  measurement  of  y^(t()) 
while  the  other  measurement  was  corrupted  by  +  em.  Point  3  depicts  the  case  where 
both  measurements  were  affected  by  the  same  positive  error.  In  the  case  of  point  4. 
'2^0*  was  measured  with  no  error  but  the  error  in  Vj(tQ)  is  +em.  Point  5  is.  in  a  way, 
symmetric  to  point  1:  yiCtn)  is  measured  corrupted  by  its  maximum  allowable  error 
while  yoltf))  ^s  measured  with  its  minimum  allowable  bound.  In  the  case  of  point  6.  no 
error  is  assumed  while  measuring  yi(tQ),  but  the  second  measurement  of  the  output  is 
affected  by  an  error  equal  to  -em.  Point  7  is  symmetric  to  point  3:  both  measurements 
of  the  output  are  affected  by  the  same  negative  error  -em.  Finally,  point  S  is  the  case 
where  no  error  was  present  in  the  measurement  of  }'->{Iq)  but  }'j(tQ)  was  affected  by  an 
error  equal  to  -e  .  In  all  cases  the  solution  of  the  system  of  equations  is  given  by  the 
projection  of  these  points  onto  the  axes  Xj(tQ)  and  x^Uq). 
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Figure  3.1     Graphical  interpretation  of  ill-conditioned  svstems-I. 
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The  discussion  above  gives  enough  information  to  answer  the  question  of  how 
the  true  solution  of  the  system  of  equations  is  perturbed  when  errors  in  the 
measurements  exist.  As  seen  in  Figure  3.1.  the  projections  of  points  1  and  5  onto  the 
axes  yield  the  most  perturbed  solutions.  On  the  other  hand,  the  projection  of  points  3 
and  7  yield  solutions  that  are  as  close  as  possible  to  the  true  solution,  assuming  that,  if 
any  errors  are  made  in  the  measurements,  they  are  equal  to  ±e  .  Therefore,  if  the 
errors  in  the  measurements  are  of  opposite  signs  and  equal  magnitudes,  the  most 
perturbed  solutions  are  obtained.  In  other  words,  if  the  errors  in  the  measurements  are 
as  far  as  they  can  possibly  be  from  one  another  the  corresponding  solution  of  the 
system  of  equations  is  also  as  far  as  it  can  be  from  the  true  solution.  Conversely,  if  the 
errors  in  the  measurements  are  equal  both  in  magnitude  and  in  sign,  the  least 
perturbed  solution  for  this  maximum  error  case  is  obtained.  Thus,  points  1  and  5 
correspond  to  the  worst  possible  case  when  the  errors  corrupting  the  measurements  are 
maximum  while  points  3  and  7  correspond  to  the  best  possible  case  when  the 
maximum  errors  in  the  measurements  are  present.  The  projection  of  points  number  1 
and  5  onto  the  axes  Xi(tn)  and  x^Ctn)  gives  the  maximum  possible  uncertainty  region 
for  the  solution  of  the  system  of  equations.  As  can  be  seen  in  Figure  3.1.  these 
uncertainty  regions  are  of  the  same  order  of  magnitude  for  both  states  but  much  larger 

than  the  errors  made  in  the  measurements. 

T 

Suppose  that  the  Q     matrix  produces  lines  that  are  perpendicular.  As  can  be 

seen  in  Figure  3.2  the  points  1,  5,  3  and  7  produce  solutions  that  are  nearly  equally 
disturbed.  Furthermore,  the  uncertainty  regions  in  the  solution  of  the  system  of 
equations  are  of  the  same  order  of  magnitude  as  the  errors  in  the  measurements. 
Therefore,  for  the  cases  where  the  lines  describing  the  system  are  orthogonal,  the 
uncertainty  regions  in  the  solution  of  the  system  of  equations  are  of  the  same  order  of 
magnitude  as  the  errors  made  while  measuring  the  output. 

Notice  the  similarities  and  differences  between  the  cases  depicted  in  Figures  3.1 
and  3.2.  In  both  cases  the  maximum  errors  in  the  measurements  are  assumed  to  be 
equal.  However,  in  the  case  of  Figure  3.1  the  errors  in  the  measurements  are 
enormously  amplified  when  the  computation  of  the  states  takes  place,  while  in  the  case 
of  Figure  3.2  the  amplification  is  practically  nonexistent.  In  both  cases  the  uncertainty 
regions  for  both  state  variables  are  nearly  equal  yielding  the  conclusion  that  both  state 
variables  are  equally  observable. 


23 


Figure  3.2     Graphical  interpretation  of  ill-conditioned  systenis-II. 
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What  is  behind  the  difference  oi"  behaviour  in  these  two  cases?  One  thing  that 
rapidly  comes  to  mind  is  the  relative  position  o(  the  lines  describing  the  system  under 
study.  In  the  case  depicted  in  Figure  3.1,  these  lines  have  almost  the  same  slope,  i.e.. 
they  are  nearly  parallel.  In  the  case  of  Figure  3.2  the  lines  are  orthogonal. 

In  the  cases  studied  so  far  both  state  variables  had  approximately  the  same 
degree  of  observability.  What  happens  if  the  states  have  different  degrees  of 
observability?  Figure  3.3  depicts  such  a  case.  Suppose  that  the  QT  matrix  produces 
lines  whose  angles  with  the  Xj(tQ)  axes  are  very  small  when  compared  with  the  angles 
between  the  lines  and  the  X2(t0)  axis.  In  this  case  the  projection  of  the  points  1  and  5 
onto  the  axes  Xj(tQ)  produces  an  uncertainty  region  that  is  much  larger  than  the  errors 
in  the  measurements.  The  projection  of  the  same  points  onto  the  ^^o)  axes  produces 
an  uncertainty  region  that  is  much  smaller  than  the  one  obtained  for  the  other  state 
variable  and  it  is  also  of  about  the  same  order  of  magnitude  as  the  error  in  the 
measurement.  Thus,  the  same  error  made  in  the  measurements  is  enlarged  a  substantial 
amount  when  the  computation  of  Xj(tQ)  takes  place  while  it  remains  unchanged  when 
X2(tQ)  is  computed.  Therefore,  according  to  the  definition  of  relative  observability  one 
can  say  that  state  X2Uq)  is  relatively  more  observable  than  state  Xi(t(j).  Notice, 
however,  that  if  by  any  chance  the  errors  in  the  measurements  are  maximum  but  with 
equal  sign  (  points  3  and  7  ),  the  uncertainty  region  for  Xj(tA)  is  very  small  while,  for 
X2(tA>,  the  uncertainty  region  doesn't  differ  much  from  the  case  of  points  1  and  5. 
Therefore,  the  sensitivity  of  the  solution  to  the  variation  of  the  errors  is  also 
meaningful  insofar  as  the  relative  observability  is  concerned:  if  the  state  is  weakly 
observable  it  is  also  more  sensitive  to  changes  in  the  errors  of  the  measurements.  If  the 
state  is  strongly  observable  it  is  not  very  sensitive  to  changes  in  the  errors  of  the 
measurements. 

The  three  cases  described  so  far  revealed  two  kinds  of  information:  one  imbedded 
in  the  relative  position  of  the  lines  describing  the  system  and  the  other  imbedded  in  the 
position  of  these  lines  with  respect  to  the  coordinate  axis.  The  first  one  determines  the 
order  of  magnitude  of  the  uncertainty  regions  of  the  computed  state  variables:  lines 
nearly  orthogonal  imply  uncertainty  regions  for  the  states  of  the  same  order  of 
magnitude  as  the  uncertainty  regions  for  the  measurements.  Lines  nearly  parallel  imply 
the  possibility  of  at  least  one  uncertainty  region  for  the  states  much  larger  than  the 
uncertainty  resion  for  the  measurements. 


Figure  3.3     Case  of  different  degrees  of  observability. 
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The  second  information  is  the  more  important  because  it  seems  to  be  directly  related  to 
the  evaluation  of  the  relative  degree  of  observability  for  each  state  variable.  In  the 
cases  depicted  in  Figures  3.1  and  3.3  the  relative  positions  of  the  lines  were  nearly  the 
same.  However,  their  positions  with  respect  to  the  coordinate  axes  were  not:  in  the 
case  depicted  in  Figure  3.1  the  angles  between  the  lines  and  both  axes  were  about  the 
same  while  in  the  case  depicted  in  Figure  3.3  the  angles  between  the  lines  and  Xj(t()) 
were  very  small  and  the  angles  between  x2(t())  were  a^most  90  degrees.  As  previously 
discussed,  in  the  case  of  Figure  3.1  both  state  variables  were  equally  weakly  observable. 
In  the  case  of  Figure  3.3  state  variable  Xi(t/\)  was  weakly  observable  while  x->(tQ>  was 
relatively  strongly  observable.  Thus,  the  angles  between  the  lines  and  the  coordinate 
axes  are  strongly  related  to  the  relative  degree  of  observability  of  the  state  variables:  if 
these  angles  are  small  (  lines  nearly  parallel  to  the  axes  ),  the  corresponding  state 
variable  is  weakly  observable.  If  the  angles  are  nearly  90  degrees  the  corresponding 
state  variable  is  strongly  observable. 

Figure  3.4  shows  the  case  depicted  in  Figure  3.1  deprived  of  several  auxiliary  lines 
to  enable  a  better  understanding  of  the  discussion  to  follow.  Consider  the  lines 
perpendicular  to  those  describing  the  system:  the  equation  of  the  line  orthogonal  to 
the  line  described  by  equation  3.20  that  passes  through  the  origin  is  given  by  : 

1.0lx2(t0)  =  0.99x1(t0)  (3.22) 

The  equation  of  the  line  orthogonal  to  the  line  described  by  equation  3.21  that 
passes  through  the  origin  is  given  by  : 

1.01x1(t0)=0.99x2(t0)  (3.23) 

Consider,  now,  the  unit  vectors  Uj  and  u2  such  that  Uj  e  x^t^)  and  u2  e  x2(t0).  As 
can  be  seen  in  Figure  3.4  ,  the  vectors  Vj  and  v-,  belong  respectively  to  the  lines 
described  by  equations  3.22  and  3.23  and  are  given  by: 

v1  =  1.01u1+0.99u2  (3.24) 

and 

v2=0.99u1+1.01u2  (3-25) 
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Figure  3.4     Relation  between  the  row  entries  of  Q 
and  the  vectors  v>  and  v^. 


T 
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The  components  of  the  vectors  Vi  and  v->  are  the  entries  in  the  first  and  second 

T 

rows  of  Q1  or  the  first  and  second  columns  of  Q.  As  can  be  seen  in  Figure  3.4  the 
angles  between  the  vectors  \j  and  \^  and  the  coordinate  axes  are  nearly  equal,  i.e., 
P  j      P2-  This  Figure  depicts  the  case  where  both  state  variables  have  equal  degrees  of 

observability.    Thus,  equal  angles  between  the  vectors,  whose  components  are  the  row 

T 

entries  of  the  Q     matrix,  and  the  coordinate  axes  may  imply  that  both  state  variables 

are  equally  observable. 

Figure  3.5  shows  the  case  depicted  in  Figure  3.3  where  state  variable  x-,(t(j)  was 
relatively  more  observable  than  state  variable  Xj(tQ).    In  this  case  both  Vi  and  v->  are 

"closer''  to  ^(tg),  LQ->  ®\  <®2  anc*  ^1  <^2~ 

Analyzing  and  comparing  the  two  cases  just  discussed,  the  importance  of  the 
information  imbedded  in  the  angles  between  the  vectors,  (  whose  components  are  the 
row  entries  of  the  Q1  matrix  ),  and  x^q)  and  X2(t0),  can  be  understood:  the  smaller 
the  angle  between  the  normal  vectors  Vj  and  Wy  anc*  an  axis  of  the  state-space,  the 
higher  the  accuracy  of  computing  the  state  corresponding  to  that  axis.  i.e.  the  higher 
the  degree  of  observability  for  that  state  variable.  The  angles  between  the  vectors  and 
the  coordinate  axes  are  directly  related  to  the  magnitudes  of  their  components  along 
these  axes:  the  larger  the  magnitude  of  the  component  along  a  particular  axes,  the 
smaller  the  angle  between  the  vector  and  the  axis.  As  previously  discussed,  the  smaller 
the  angle  the  higher  the  accuracy  when  computing  that  state.  Recalling  that  the  rows 
of  Q  are  the  vectors  relating  Xj(tQ)  and  x->(tQ)  with  VjUq)  and  y->(tQ).  an  easy  way  to 
check  the  relative  observability  of  a  state  variable  is  to  look  for  the  largest  principal 
axis  component  of  each  normal  vector  and  then  associate  these  components  with  the 
state  variables.  This  idea  has  been  used  as  an  algorithm  for  the  evaluation  of  the 
relative  degree  of  observability  of  a  system  and  is  one  of  the  three  algorithms  presented' 
in  the  Fortran  program  TOBS  in  Appendix  A. 

F.       A  MORE  GENERALIZED  METHOD 

The  graphical  method  previously  described  is  quite  helpful  for  low-order  systems 
because  it  requires  few  operations,  is  rapid  and  gives  a  good  feel  for  the  rank  of  the 
state  variables  as  far  as  relative  observability  is  concerned.  When  the  system  becomes 
larger  and  more  complex  the  direction  of  minimum  error,  although  perpendicular  to  the 
direction  of  maximum  error,  may  be  difficult  to  find  due  to  the  increasing  number  of 
vectors  in  a  higher  order  state-space.   What  if  instead  of  checking  the  whole  space 


Fieure  3.5     Importance  of  the  anelcs  on  evaluating 
"the  degree  of  observability  per"state- variable. 
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T 
spanned  by  Q     one  finds  a  way  of  comparing  one  particular  vector  against  the  state- 
space  axes  to  obtain  the  relative  degree  of  observability  for  each  state  variable  ? 

The  direction  of  minimum  error  conveys  considerable  information  about  the 
relative  observability  of  the  system.  Let  W|.w->....w  be  the  columns  of  Q  normalized. 
In  [Ref  7:  p.  9]  the  author  forms  an  observability  function  called  \\i  which  is  a  scalar 
that  is  given  by 

\|/  =  (w1Tu)2  +  (w2Tu)2  +  ...  +  (wnTu)2  (3.26) 

where  u  represents  a  vector  in  the  direction  of  minimum  error.  The  length  of  this 
vector  is  constrained  to  be  of  unit  length,  i.e., 

uTu=l  (3.27) 

Equation  3.26  can  be  written  as  follows  : 

y  =  uT[w1w1T  +  w2w2T..  +  wnwnT]u  (3.28) 

The  quantity  inside  the  brackets  in  equation  3.28  is  actually  QnQn  [Ref.  12:  p.  220]. 
With  this  result  in  mind,  this  last  equation  can  be  written  as 

y  =  uT(QnQnT>u  (3.29) 

The  goal  o[  this  derivation  is  to  find  a  vector  that  gives  the  direction  oi" 
maximum  error,  yielding  a  minimum  for  the  observability  function  i|/  .  To  solve  this 
problem,  therefore,  the  gradient  of  vj/  with  respect  to  u  must  be  found.  In  this  problem 

there  is  a  constraint  that  is  given  by  equation  3.27.    In  general,  any  quadratic  form  of 

T 

the  type  x   Px  (  P  positive  definite  )  is  bounded  below  by  the  product  of  the  minimum 

eigenvalue  of  P  and  the  length  of  x  and  bounded  above  by  the  product  of  the 
maximum  eigenvalue  of  P  and  the  length  of  x  [Ref.  13:  p.  121].  For  the  particular  case 
of  equation  3.29,  recalling  that  the  length  of  u  was  assumed  to  be  unity,  this  result  can 
be  written  as 

Xrmn^uT(QnQnT>u^max  <3'30) 
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where 

^min  *s  ^  sma^esl  eigenvalue  of  QnQn 

T 
^max  *s  t^ie  ^ar§est  eigenvalue  of  QnQn 

The  same  result  can  be  obtained  using  direct  differentiation  of  equation  3.29  with 
respect  to  u  [Ref.  7:  p.  10].  In  this  procedure  the  Lagrangian  multiplier  formulation  is 
used  where  the  constraint  is  given  by  equation  3.27.  The  Lagrangian  multiplier 
formulation  is  given  by  equation  3.31. 

— [uT(QnQnT)u  -  X(uTu-l)]  =  0  (3.31) 

du 

The  details  on  the  differentiation  of  equation  3.31  can  be  found  in  [Ref.  14:  p.  288]. 
The  result  of  the  differentiation  is  given  by  equation  3.32 

(QnQnT-Xl)u  =  0  (3.32) 

T 

Rearranging  equation   3.32   and  premultiplying   both   sides   by  u      equation   3.33   is 

obtained. 


u1(QnQnI)u=uIXu=X  =v|/  (3.33) 

This  equation  shows  that  the  minimum  value  of  the  observability  function  \\f  is  equal 
to  the  smallest  eigenvalue  of  QnQn  •  Therefore,  the  direction  of  maximum  error  ( 
where  the  observability  function  is  minimum  )  will  be  given  by  the  eigenvector 
associated  with   the   smallest  eigenvalue.     The  algorithm  for  the  evaluation   of  the 

relative  observability  for  each  state  variable  is,  now,  evident  and  is  as  follows: 

T 

1.  form  Qn  ,  Qn     and  multiply  them  together  -  call  the  product  matrix  QQT 

2.  find  the  eigenvalues  and  associated  eigenvectors  of  QQT 

3.  find  the  smallest  eigenvalue  and  associated  eigenvector 

4.  sort  the  components  of  the  eigenvector  associated  with  the  smallest  eigenvalue 
in  increasing  order  of  magnitudes 

5.  associate  the  sorted  components  with  the  respective  state  variables 

6.  the  smallest  component  is  associated  with  the  most  observable  state  while  the 
largest  component  is  associated  with  the  least  observable  state 
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G.      SUMMARY 

In  this  Chapter  the  concept  of  relative  degree  of  observability  was  given  as  a 
measure  of  the  accuracy  involved  in  determining  the  initial  states  from  the 
measurements  of  the  output.  The  concept  of  an  ill-conditioned  matrix  was  presented  as 
well  as  the  problem  involved  in  inverting  highly  ill-conditioned  matrices.  Based  on  the 
definition  of  relative  degree  of  observability  given  and  the  description  of  the  problem  of 
an  ill-conditioned  matrix,  four  methods  for  evaluating  the  relative  degree  of 
observability  were  presented:  the  upper-bound  error  method,  the  standard  deviation 
error  method,  the  graphical  method  and  the  generalized  graphical  method.  Examples 
were  presented  for  the  first  two  methods  and  algorithms  were  derived  for  the  first,  third 
and  fourth  methods.  These  algorithms  are  all  imbedded  in  the  same  Fortran  program 
(TOBS)  presented  in  Appendix  A. 
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IV.  IMPROVEMENT  OF  THE  DEGREE  OF  OBSERVABILITY  OF  A 
WEAKLY  OBSERVABLE  SYSTEM 

A.  INTRODUCTION 

In  the  last  Chapter,  methods  for  evaluating  the  degree  of  observability  for  each 
state  were  presented.  A  natural  extension  of  the  study  of  this  problem  is  the  search  for 
a  method  of  improving  the  observability  of  a  particular  weakly  observable  state.  In  this 
Chapter  methods  for  improvement  will  be  presented,  both  for  the  single  and  multi- 
output  cases. 

B.  THEORETICAL  BACKGROUND 

Before  attempting  to  derive  the  method  to  increase  the  relative  observability  of  a 
state  variable,  some  important  concepts  of  linear  algebra  need  to  be  reviewed. 
Consider,  then,  the  system  represented  by  equation  2.1  where  the  square  matrix  A  is 
such  that  all  its  eigenvalues  are  distinct.  In  this  case  the  matrix  A  has  at  least  one 
eigenvector  for  each  distinct  eigenvalue  that  satisfies  equation  -4.1 

Ap=Xp  (4.1) 

where  p*0  is  the  eigenvector  of  the  square  matrix  A  associated  with  the  eigenvalue  X 
{X  is  a  scalar)  [Ref.  10:  p.  661].  Since  the  assumption  is  that  the  eigenvalues  of  A  are 
all  distinct,   the  n  corresponding  eigenvectors  will  all  be   linearly  independent.    Let 

P=[pj p  ]   be   the  matrix   of  the  eigenvectors   of  A.     The   matrix   M  =  P      exists 

because  the  columns  of  P  are  linearly  independent.  Let  m^     be  the  rows  of  M.  It  can 

T  T  T  T 

be  shown  that  A  1rr^=X-m^  or  equivalently  m^  A=>,1m1  .  The  m^  are  often  called 
the  left  (  or  row  )  eigenvectors  of  A  in  contrast  to  the  right  (  or  column  )  eigenvectors 
of  A  [Ref.  10:  p.  664].    Note,  however,  that  the  so  called  left  eigenvectors  of  A  are 

T 

nothing  but  the  regular  eigenvectors  of  A  .  In  [Ref.  10:  p.  664]  the  author  shows  two 
relations  that  will  be  quite  important  in  the  derivation  of  the  method  described  in  the 
next  section.  These  relations  are  described  as  being  the  spectral  decomposition  of  A 
and  are  siven  bv 
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=  PAM 
=  M"1AM 


(4.2) 


\2     v  \2         T 

A  =L  x  iPi1^ 

=  PA2M 

=  M'1A2M 


(4.3) 


In  general  the  spectral  decomposition  of  An  is  given  by 


An  =  IXnipimi 
=  PAnM 

=  \rIAnM 


I 


(4.4) 


where  A=  diag  [X.j,X2,.-M^n}-   The  relations  presented  in  the  equations  4.2,  4.3  and  4.4 
will  be  very  helpful  in  understanding  the  method  that  is  presented  in  the  next  section. 

C.       DERIVATION  OF  THE  METHOD 

Consider,  now.  the  transpose  of  the  observability  matrix  defined  in  equation  2.3  : 


qt= 


C 
CA 

CA2 


CA 


n-1 


(4.5) 


Substituting  the  representation  of  the  matrix  A  given  by  equations  4.2  and  4.3,  yields 
equation  4.6  : 


QT  = 


CM_1M 

CM_1AM 
CM_1A2M 


CM^A^M 


(4.6) 
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Notice  that  in  equation  4.6  the  matrices  C.  M      and  M  are  common  to  all  the  rows. 
Letting 


L  =  CM"1, 


(4.7) 


substituting  equation  4.7  in  equation  4.6  and  post-factoring  the  matrix  M,  equation  4.8 

is  obtained  : 


QT  = 


L 
LA 

LA: 


LA 


n-1 


M 


(4.8) 


The  matrix  L,  is  a  row  vector  of  dimension  lxn,  because  C  is  a  lxn  row  vector  (  in  the 
case  of  single  output  )  and  M      is  a  nxn  matrix.  Thus,  L  can  be  written  as 


L=[L1L2...Ln] 


(4.9) 


With  this  in  mind  and  recalling  that  that  A  is  a  diagonal  matrix  whose  elements  are  the 

T 

eigenvalues  of  A.  the  matrix  Q     can  be  rewritten  as  follows  : 


0 


T_ 


Ll>-12 


L2X2 
L2X2- 


Lpi  ^2   2 


Ln*x 
LnV 


LnK 


n-1 


M 


(4.10) 


,T 


Looking  carefully  at  equation  4.10  it  can  be  seen  that  the  matrix  Q     can  be  written  as 
the  product  of  three  matrices  : 


QT=VGM 


(4.1L 
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where  V  is  the  Vandermonde  matrix  of  the  eigenvalues  of  A: 


V    = 


1 


X  n-1 


1 

^2 


\    n-1 


n-1 


(4.121 


G  is  a  diagonal  matrix  whose  elements  are  the  entries  of  L: 


G  = 


0 

u 


0 
0 


(4.13; 


and  M  is  the  matrix  of  the  left  eigenvectors  of  A. 

T 

With  this  procedure  the  matrix  Q     has  been  transformed  into  a  product  of  three 

matrices.  Since  V  and  M  depend  only  on  the  plant  and  G  depends  on  both  the  plant 
and  the  C  matrix,  it  seems  to  be  G  that  needs  to  be  studied  and  changed  in  order  to 

improve  the  relative  degree  of  observability  of  a  particular  state  variable. 

T 

It   was   mentioned   in   Chapter   III    that   the   degree   of  ill-condition   of  Q      is 

important  in  determining  the  magnitude  of  the  errors  (ec)  in  computing  the  initial 

states  once  an  error  is  introduced  in  measuring  the  output  (em).    Therefore,  it  is 

T 
valuable  to  have  a  quantitative  measure  of  the  degree  of  ill-condition  of  the  matrix  Q   . 

In  [Ref.  15:  p.   121]  the  author  defines  the  condition  number  of  a  matrix  W  as  the 

product  of  two  matrix  norms  : 


cond(W)=||W||  ||\V 


rli 


(4.14) 


If  the  condition  number  is  small  the  matrix  W  is  not  ill-conditioned.  If  the  matrix  W  is 
ill-conditioned  the  entries  of  its  inverse  will  be  large  numbers,  yielding  a  large  condition 
number  for  W.  However,  this  can  also  be  true  when  the  elements  of  W  are  small  even 
in  the  absence  of  ill-conditioning.  Multiplying  the  two  norms  has  a  normalization 
effect,  so  the  condition  number  is  large  only  for  an  ill  conditioned  system  [Ref.  15:  p. 
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122].  Observe  that  the  condition  number  cannot  be  smaller  than  1,  which  corresponds 
to  the  condition  number  of  the  identity  matrix  [Ref.  15:  p.  121]. 

Equations  3.3  and  3.4  show  that  the  vector  of  the  measurements  of  the  output 
consists  of  the  sum  of  the  actual  values  of  the  output  [ya(tn)]  with  the  error  (e  ) 
introduced  while  measuring  the  output.  In  other  words  y(i(tQ)  =  ya(tn)  +  e  ,  which  is 
equivalent  to 

em=>'d(t0)-ya(t0)  <4-15) 

T  T 

But    ya(t())  =  Q   xa(t())    ant*   y(j(to)=Q   x(to)-     By    substituting    these    expressions    in 

equation  4.15.  equation  4.16  is  obtained  : 

em=QTx(t0)-QTxa(t0)  (4.16) 

=  QT[x(t0)-xa(t0)] 
=  QTec 

The  error  in  the  computed  states  is  given  by 

ec  =  (QT)"lem  <4-17) 

Taking  norms  and  recalling  that,  in  general,  ||AB||  ^  ||A||  ||B||,  yields  for  the  particular 
case  of  equation  4.17  the  result 

l|ecll  <  IKQ"1-)-1!!  ||em||  (4.18) 

Applying  norms  to  the  last  part  of  equation  4.16  the  following  is  obtained  : 

HemII^IIQT||  ||ec||  (4.19) 

Notice  that  equations  4.18  and  4.19  constitute,  respectively,  upper  and  lower  bounds 
on  the  magnitude  of  the  error  in  computing  the  initial  states.  This  fact  is  shown  in 
equation  4.20 
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^l|ec||^||(QT)-1||  ||em||  (4.2.), 


Q 


Applying  the  same  reasoning  to  the  equations  relating  xa(tQ>  to  yjtn),  yields  for 
xa(tQ>  the  upper  and  lower  bounds  shown  in  equation  4.21, 

iiya<to)ii  t  1 

-j^-^IIXa^ll^lKQ7)-1!!!^^)!!  (4.21) 

which  is  equivalent  to  the  result  obtained  in  equation  4.22 

l  l  IIQTII 

— rn < <  (4.22) 

ll(QT)  l\\  llya(t0)ll       Hxa(t0)||      ||ya(t0)|| 

Multiplying  this  last  inequality  by  the  inequality  4.20  and  recalling  the  definition  of 
condition  number  given  in  equation  4.14.  equation  4.23  is  obtained  for  the  bounds  of 
the  relative  error  in  computing  the  initial  states. 


1  lie    II  lie  II  lie    II 

-  <  (cond  QT)        m  (4.23) 


(condQ1)     ||ya(t())||      ||xa(t0)|]  ||ya(t0) 

This  last  result  shows  that  the  relative  error  in  the  computed  solution  (eJ  can  be  as 
large  as  the  relative  error  in  the  measurements  multiplied  by  the  condition  number.  Of 
course  it  can  also  be  as  small  as  the  relative  error  in  the  measurements  divided  by  the 
condition  number.  So,  when  the  condition  number  is  large,  the  error  in  the 
measurements  gives  little  information  about  the  accuracy  of  x(tQ).  Conversely,  if  the 
condition  number  is  nearly  unity,  the  relative  error  in  the  measurements  is  a  good 
measure  of  the  relative  error  in  computing  the  states  [Ref.  15:  p.  123],  and  the  error  in 
the  measurements  and  the  error  in  the  calculated  states  are  of  the  same  order  of 
magnitude.  The  goal  of  this  whole  procedure  is  to  improve  the  relative  degree  of 
observability  of  a  weakly  observable  state.  Therefore,  if  the  condition  number  of  Q  is 
minimized  the  relative  error  in  the  measurements  will  be  of  the  same  order  of 
magnitude  as  the  error  in  computing  the  states.  The  optimum  result  will  be  obtained 

for  a  condition  number  equal  to  one.  In  this  case  the  upper  and  lower  bounds  on  ec 

T 

will  coincide.  The  condition  number  of  Q     is  given  by 
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cond(QT)=||QT|i  IKQ1)"1!!  (4.24) 

T 

Recalling  equation  4.11  where  Q    =  VGM  and  substituting  this  in  equation  4.24  yields 

cond(QT)  =||VGM||  ||(VGV[)_1||  (4.25) 

=  MVGMII  IIM^G'H'-1]! 


But  the  norm  of  the  product  of  three  matrices  must  not  exceed  the  product  of  the  three 
norms,  i.e., 

cond(QT)  <||V||  ||G||  ||M||  ||\rl||  NG"1}!  i|Vl||  (4.26) 

^(cond  V)(cond  G)(cond  VI) 

Observe  that  V,  V  ,  M  and  M  cannot  be  changed  because  they  are  totally 
dependent  on  the  dynamic  characteristics  of  the  plant.  Therefore,  only  by  minimizing 
the  condition  number  of  G  will  the  condition  number  of  Q     be  minimized.    The  matrix 

G  is  a  function  of  the  entries  of  C  (equations  4.7  and  4.13),  which  are  the  only  parts  of 

T 

the  system  that  the  user  can  change.  On  the  other  hand  the  condition  number  of  Q     is 

a  function  of  the  condition  number  of  G  (equation  4.26).  Therefore,  the  relation  that 
exists  between  the  condition  number  of  G  and  the  condition  number  of  Q     can  be  seen 

as    a  function  whose  independent  variable  and  abcissa  is  the  condition  number  of  G 

T 

and  the  ordinate  is  the  condition  number  of  Q   . 

(condQT)=f(condG)  (4.27) 

The  goal  of  this  procedure  is  to  find  the  minimum  of  the  condition  number  of  Q  by 
changing  a  specific  entry  in  the  C  matrix.  In  other  words,  the  goal  is  to  find  the 
minimum  of  f.  In  order  to  obtain  an  algorithm  that  accomplishes  this,  the  function  f 
must  be  studied  in  more  detail.  Both  the  abcissa  and  the  ordinate  are  condition 
numbers  of  matrices  which,  by  definition,  can  never  be  less  than  one.  Only  in  very 

special  cases  is  the  condition  number  of  a  matrix  equal  to  one.  The  identity  matrix  and 

T 

all  its  multiples  are  among  those  special  cases.  Notice,  also  that  the  first  row  of  Q     is 

the  C  matrix  in  the  single  output  case.  Therefore,  only  when  c  1 1  is  different  from  zero 
and  all  the  other  entries  in  the  C  matrix  are  zero,  will  the  possibility  of  obtaining  a 
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minimum  for  the  condition  number  of  Q  equal  to  one  be  likely  to  occur.  In  other 
words,  in  the  majority  of  cases  both  abcissa  and  ordinate  will  be  greater  than  one.  The 
matrix  G  plays  an  important  role  in  the  development  of  the  algorithm  that  minimizes 
the  condition  number  of  Q1  by  changing  a  specific  entry  in  the  C  matrix.  In  order  to 
understand  how  the  matrix  G  is  used  in  this  derivation  a  more  detailed  study  of  its 
entries  will  be  done.   The  matrix  G  is  a  diagonal  matrix  whose  elements  are  given  by 

Ll  =  §ll=Iclj^l 
L2  =  §22  =  Icljnii2 

Ls=gSs  =  Icljmis 

•      .      .      .  j— 1,2 n  (4.2S) 

Lf  =gff  =Icljm]€ 

Ln=§nn  =  Icljny 

where  the  Cj:  are  the  entries  of  the  C  matrix  and  the  m:  are  the  entries  of  the  M 
matrix. 

The  norm  used  to  compute  the  condition  number  may  be  any  matrix  norm 
[Ref.  15:  p.  118].  In  the  present  case  the  maximum  column  sum  norm  will  be  used  to 
compute  the  condition  number  and  will  be  defined  as  follows  : 

IIGH^maxYlgjjl         j=  1,2 n  (4.29) 

Notice  that  both  subscripts  of  g  are  equal  because  the  matrix  G  is  diagonal.  In  this 
particular  case  the  sum  is  equal  to  the  only  non-zero  entry  in  the  column.  The  norm  is 
then  given  by  equation  4.30: 

[|G||1  =  max|gjj|  j=  1,2 n  (4.30) 

Suppose  that  the  maximum  element  in  G  is  gf£  where  g^ p  is  given  by  equation  4.28. 
According  to  equation  4.30  the  norm  of  G  will  then  be  equal  to 

MOM!  =  |g££|  (4.31) 
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Since  the  G  matrix  is  diagonal  its  minimum  entry  corresponds  to  the  maximum  entry 
in  G     which  is  also  diagonal.  Its  norm  is  given  by  equation  4.32  : 


or 


IG'1!!-       *  i— 1.2 n  (4.32) 

minlgjjl 


!G-'ll=  (4.33) 


where  gss  is  the  smallest  entry  in  G  and  is  given  by  equation  4.28  .  The  condition 
number  of  G  is,  by  definition,  given  by  the  product  of  the  two  norms  given  in 
equations  4.31  and  4.33.  Substituting  equation  4.28  for  the  values  of  gss  and  g££  in 
equations  4.31  and  4.33  yields,  for  the  equation  of  the  condition  number  of  G 

n 


n 

S 


■=  (condG)  (4.34) 


3  =  1 


The  ratio  given  by  equation  4.34  is  also  the  ratio  between  the  largest  and  smallest 
eigenvalues  of  G  (G  is  a  diagonal  matrix),  which  is  another  definition  of  condition 
number  of  a  matrix  [Ref.  16:  p.  38].  The  right-hand  side  of  equation  4.34  is  a  constant. 
Let  the  inverse  of  that  constant  be  equal  to  K.  Equation  4.34  can  then  be  rewritten  as 
follows: 

n  n 

j=l  j=l 
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Since  there  is  only  one  equation  (equation  4.35  )  to  provide  the  changes  in  the  C 
matrix  to  improve  the  relative  degree  of  observability  of  the  system,  a  decision  has  to 
be  made  upon  the  entry  in  the  C  matrix  for  which  equation  4.35  will  be  solved. 
Assuming  that  at  this  stage  one  has  already  computed  the  relative  degree  of 
observability  of  the  system,  one  knows  which  is  the  least  observable  state.  If  that  state 
is  state  xr  equation  4.35  will  be  solved  for  Cjr  because  this  is  the  entry  in  the  C  matrix 
associated  with  the  least  observable  state.  Expanding  equation  4.35  yields  : 

K(cnml£  4-...  +  clrmr€  +  ...  +  clnmn|  )=  (4.36) 

c11mls  +  ...4-cIrmrs+...  +  clnmns 

Rearranging  and  collecting  common  terms. 

(Kmu-mls)c11  +  ...  +  (Kmrrrnrs)clr+...  +  (Kmnrnins)cln=  0  (4.37) 

Solving  for  Cjr  equation  4.3S  is  obtained 

clr=      -4-lPiclr         Pr*°  V-™ 

^r    i=l 
i*r 

where  in  general  pj  =  (Kmj£  -m:.).    At  this  point  every  tool  is  available  to  develop  the 

algorithm  that  will  provide  the   new  entry   in  the   C   matrix  corresponding  to   the 

T 

minimum  condition  number  of  Q   . 

The  initialization  of  the  algorithm  is  done  by  selecting  an  abcissa  equal  to  the 

condition  number  of  G  of  the  initial  system  and  an  ordinate  equal  to  the  condition 

T 

number  of  Q     of  the  same  system.    To  define  the  interval  where  the  minimum  is  to  be 

evaluated,  another  point  has  to  be  chosen.  As  previously  mentioned,  the  condition 
number  of  a  matrix  is  a  quantity  that  is  always  greater  or  equal  to  one.  With  this  in 
mind,  the  choice  of  the  second  point  selected  for  the  initialization  of  the  algorithm  is 

the  one  where  the  abcissa  (condition  number  of  G)  is  equal  to  one  and  the  ordinate  is 

T 

the  corresponding  value  for  the  condition  number  of  Q    .  The  computation  of  the 

T 

condition  number  of  Q     given  the  condition  number  of  G  is  a  three  step  procedure 

that  can  be  summarized  as  follows: 

STEP-1  :  Given  the  condition  number  of  G  (abcissa)  solve  equation  4.38  for  Cjr 
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STEP-2  :  I-ind  the  new  Q1  by  using  equation  2.3  and  transposing  it. 

T 
STEP-3  :  Find  the  condition  number  of  the  new  Q    (new  ordinate). 

A  flow-chart  of  the  algorithm  with  the  initialization  and  minimization  procedure 

is  presented  in  Figure  4.1.  In  the  flow-chart  the  letter  f  corresponds  to  the  definition  in 

equation  4.27  and  represents  the  three  step  procedure  described  above.  X  represents 

the  condition  number  of  G  while  Y  represents  the  condition  number  of  Q    .    After 

bisecting  the  interval  a  decision  has  to  be  made  concerning  where  to  search  for  the 

minimum  of  f.  i.e..  if  to  the  right  or  to  the  left  of  X3.  The  information  given  by  the 

derivative  of  f  near  X3.Y3  is  used  to  make  this  decision:  if  the  derivative  is  positive  the 

interval  to  be  searched  is  to  the  left  of  X3.  If  the  derivative  is  negative  the  interval  to 

be  searched  is  to  the  right  of  X3.  Once  the  direction  of  search  is  known,  the  points 

limiting  the  search  interval  are  renamed:  the  point  with  the  largest  ordinate  is  X1.Y1 

and  the  point  with  the  smallest  ordinate  is  X2.Y2.  If  the  stopping  criterion  is  met  the 

algorithm  stops.  Otherwise,  the  bisecting  procedure  is  executed  again.  This  algorithm 

gives  a  local  minimum.  Other  initializations  have  to   be  used  to  try  to  find  other 

minimums.   The  whole  algorithm  containing  the  required  minimization  of  the  condition 

T 
number  of  Q    ,  was  implemented  in  a  Fortran  computer  program  called  COBS  and  is 

presented  in  the  Appendix  B. 

The  choice  of  the  entry  of  the  C  matrix  associated  with  the  least  observable  state 

as  the  one  to  be  changed  is  made  here.  However,  this  issue  is  controversial  because  real 

cases  may  happen  where  the  user  cannot  change  the  C  matrix  according  to  the  least 

observable  state.    Furthermore,  an  independent  change  in  an  entry  of  the  C  matrix  may 

not  be  possible,  i.e..  a  change  in  one  entry  of  C  may  only  be  possible  if  another  entry  is 

changed.  Several  combinations  are  possible.  Therefore,  according  to  a  specific  case,  the 

user  may  have  to  try  several  allowable  changes  in  the  C  matrix  to  determine  the 

minimum  condition  number  of  Q   .  The  program  that  is  in  Appendix  B  has  the  default 

feature  of  changing  the  entry*  of  the  C  matrix  associated  with  the  least  observable  state. 

However,  any  other  entry  in  the  C  matrix  can  be  changed  in  order  to  increase  the 

degree  of  observability  of  a  specific  state  variable.  This  is  accomplished  by  following 

the  documentation  contained  in  the  listing  of  COBS. 

T 
Once  the  condition  number  of  Q     has  been  minimized  and  the  new  entries  of  the 

C  matrix  found,  the  system's  observability  has  to  be  checked.    In  equation  4.11  it  was 

T 

seen    that    Q       could    be    decomposed    into    the    product    of   three    matrices.    The 
Vandermonde  and  the  matrix  of  the  left  eigenvectors  are  full  rank  matrices  because  by 
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assumption  the  plant  has  distinct  eigenvalues.  Therefore,  if  G  is  a  full  rank  matrix  QT 
will  also  be  a  full  rank  matrix  and  the  system  will  be  completely  observable.  For  G  to 
be  full  rank  all  its  entries  must  be  different  from  zero  because  G  is  a  diagonal  matrix. 

D.       THE  MULTI-OUTPUT  CASE 

What  has  been  said  so  far  is  only  valid  for  the  single  output  case.  Since  this  is 
very  seldom  the  case,  the  multi-output  situation  will  be  considered  in  this  section. 

Consider,  again,  the  linear  time  invariant  system  described  by  equations  2.1  and 
2.2  where  C  is  a  mxn  matrix  and  m  ^  1.  The  observability  matrix  Q  is,  for  the  general 
case  of  multi-output,  given  by 

Q  =  [CT     ATCT    (AT)2CT     ...   (AT)n-eCT]  (4.39) 

where  I  is  the  rank  of  C  [Ref.  17:  p.  503].  This  new  observability  matrix  is  nxp.  where 
p  =  (n-E+  l)m.  Notice  that  if  m=  1  and  1=  1  (single  output  case),  the  matrix  Q  will  be 
square.  The  fact  that  the  matrix  Q  is  no  longer  square,  transforms  a  relatively  simple 
problem  into  a  much  more  elaborate  one.  Recalling  equation  3.2  where  the  relation 
between  the  derivatives  of  the  output  and  the  states  is  given,  if  one  attempts  to  solve 
that  system  of  equations  for  x(tQ)  when  Q  is  not  a  square  matrix,  the  concept  of 
generalized  or  pseudo-inverse  must  be  introduced.  Premultiplying  equation  3.2  by  Q 
yields 

Qyd(t0)=QQTx(t0)  (4.40) 

Solving  equation  4.40  for  x(tQ), 

x(t0)  =  (QQT)"1Qyd(to)  (4-41) 

The  matrix  (QQ"*Y*Q  is  called  the  pseudo-inverse  of  Q  [Ref.  17:  p.  57"].  The 
solution  of  equation  3.2  using  the  pseudo-inverse  is  an  approximate  solution  in  the 
least-squares  sense  [Ref.  17:  p.  578].  Let  QQ  -  R.  The  problem  is,  now.  finding  how 
invertible  the  the  R  matrix  is.  i.e.,  how  ill-conditioned  the  matrix  R  is.  One  possible 
procedure  is  the  one  used  for  the  single  output  case.  Thus,  performing  the  spectral 
factorization  of  the  matrix  A  in  the  transpose  of  equation  4.39  by  substituting  its 
representation  given  by  equations  4.2.  4.3  and  4.4  yields  equation  4.42: 
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initialize  by  selecting 

XI  =  cond  G 

Vl  =  cond  Q1 
X2=  1 
Y2=f{l) 


bissect  interval  [XI, X2] 
->  X3  -►  Y3  =  f(X3) 


compute  derivative  near  X3 


redefine  the  search  interval 

using  the  information  from  the 

derivative  such  that  the  largest 

ordinate  is  X1.Y1  and  the 

smallest  ordinate  is  X2,Y2 


Figure  4.1     Algorithm  to  provide  the  minimum  of  the  condition  number  of  Q 

by  changing  a  specific  entry  in  C. 
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Qt= 


c\rl\i 

CM"1  AM 
CM_1A2M 

CM-]An-C 


M 


(4.42) 


,-1 


Notice,  again,  that  as  in  the  single  output  case,  the  matrices  C,  M      and  M   are 
common  to  all  the  rows.  Let 


L=CM 


1 


(4.43) 


where  the  dimension  of  L  is  mxn  and 


and 


C    = 


cll 

ci2 

c21 

c22 

c31 

c32 

'ml 


'm2 


cln 
c2n 
c3n 


'mn 


(4.44) 


M'1  = 


mll 
m21 
m31 


m12 

m22 
l3 


m,2 


"ml  ^2 


m 


In 


m2n 
m3n 

"mn 


(4.45) 


Substituting  equation  4.43  in  equation  4.42  and  post-factoring  the  matrix  M.  equation 
4.46  is  obtained: 
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qt= 


L 

LA 

LA2 

LA«-£ 

M 


(4.46) 


Performing  all  the  matrix  multiplications  in  equation  4.46  except  those  involving  the 

T 

matrix  M,  the  Following  representation  for  Q     is  obtained: 


QT= 


I     i   n-£    r     \    n-6 

r     i   n-t     r     \    n-t 
'21    1  '22A2 


^lnAn 

^2nXn 


n-C 


Wri 


^m2A2 


n-C 


ua 


n-e 


mn  n 


M 


(4.47) 


where  £•:  are  the  entries  of  L  and  X>,  the  system  eigenvalues,  are  the  entries  of  the 
diagonal  of  A.  The  goal  is  to  decompose  the  first  matrix  on  the  right  hand  side  of 
equation  4.47  as  the  product  of  two  matrices  such  that  only  one  of  them  depends  on 
the  C  matrix.  The  matrix  Q1  can  actually  be  written  as  the  product  of  three  matrices: 


QT  =  VGM 


(4.48) 


where  V  has  a  Vandermonde  structure, 
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A 


A, 


A, 


V  = 


Al 


Ay 


(4.49) 


n-C 


A^ 


n-e 


A 


n-C 


Each  matrix  forming  V  is  a  mxm  matrix  of  the  form 


A>!- 


0 


0 


0 


I) 


h 


1-1 


J    J 


(4.50) 


where  i.  the  V  super-row  index,  is  equal  to  2,3,. ...n  and  j,  the  V  super-column  index  is 
equal  to  1.2,...,n.  Therefore,  V  has  in  each  of  its  super-columns  (n-£  +  1)  matrices  with 
m  rows  and  m  columns.  On  the  other  hand,  there  are  n  submatrices  in  each  super-row 
of  V.  Thus,  the  matrix  V  has  p  rows  and  q  columns  where  p  =  (n-l+  l)m  and  q  =  mn. 
The  matrix  I  is  an  mxm  identity  matrix.   The  matrix  G  is  a  qxn  block,  diagonal  matrix 


G  = 


Zq  0  0 

0  L2  0 

0  0  Ln 


(4.5V, 


Each  matrix  contained  in  G  is  actually  an  mxl  column  vector  corresponding  to  the 
columns  of  the  matrix  L  given  by  equation  4.43.  In  general,  each  submatrix  of  G  is 
siven  bv 
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k= 


<>2i 

<bi 


i=  l,2,...,n 


(4.52) 


!mi 


T 

This  decomposition  of  the  Q     matrix  seems  to  be  complicated  and  even  difficult  to  use. 

Two  examples  will  be  presented  to  clarify  the  ideas.  Suppose  that  a  third  order  plant  is 

given  in  which  only  one  output  is  available,  i.e.,  n=  3,  m=  1  and  the  rank  of  C  is  £=  1. 

T 

The  dimension  of  Q     is  pxn  where  p  =  (n-6+  l)m=  3  and  n=  3;  V  is  a  pxq  matrix  where 

T 

p=  3  and  q  =  nxm=  3;  G  is  qxn  where  q=  3  and  p=  3.  In  summary  Q   ,  V  and  G  are  3x3 

matrices.  The  matrix  V  has  in  each  of  its  columns  n-£+-  1  =  3  matrices  with  m=  1  row 
each.  According  to  equation  4.49  the  first  row  matrices  are1  identity  matrices  which 
have  m  rows  (in  this  case  one  row),  i.e.,  the  first  row  elements  are  three  ones.  The 
second  row  elements  are  found  by  using  equation  4.50  along  with  equation  4.49.  Since 
one  is  looking  for  the  second  super-row  of  V,  the  index  i  is  equal  to  2  and  the  super- 
column  index  j  is  equal  to  1.2  and  3.  Therefore,  the  matrices  that  constitute  the  second 
super-row  of  V  have  m=  1  row  and  have  entries  Xj,  \j  and  X3.  The  third  and  last 
super-row  matrices  of  V  are  Xj\  X2    and  X3". 


V    = 


1 


1 
^2 


(4.53) 


The  matrix  G  is  a  qxn,  i.e.,  3x3  diagonal  block  matrix  whose  entries  are  given  by 
equation  4.52.  Since  m=  1  the  L-  are  scalars  and  are  equal  to  L\~^\\>  L.i  =  (^\i> 
L->  =  Cj3.  The  matrix  G  will  be  equal  to 


G    = 


Cll 

0 

0 

0 

«12 

0 

0 

0 

5l3 

(4.54) 


As  can  be  seen  by  equations  4.53  and  4.54,  by  applying  the  spectral  decomposition  of 
Q  for  this  generalized  case,  a  result  is  obtained  that  agrees  with  the  one  that  would  be 
obtained  if  equations  4.12  and  4.13  were  used. 
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Suppose,  now  that  the  system  is  third  order  with  two  outputs.  In  other  words. 
n=3  and  m=2.  The  rank  of  C  is  the  smaller  of  m  and  n.  i.e..  £=2.  Thus  Q  will  be  a 
( n-C-f-  l)mxn  matrix,  i.e.,  Q  will  be  a  4x3  matrix.  This  matrix  will  be  decomposed  in 
the  product  of  three  matrices  V,  G  and  M  whose  dimensions  are  as  follows:  V  is  a  pxq 
matrix  where  p=(n-£  +  l)m  =  4  and  q  =  mxn=6.  The  matrix  G  is  qxn=6x3  while  the 
matrix  M  is  nxn=3x3.  Each  submatrix  of  V  is  mxm=2x2.  Therefore.  V  has  two  super 
rows  and  three  super  columns  because  that  is  the  only  way  of  obtaining  a  dimension  of 
4x6  for  the  matrix  V.  The  submatrices  in  the  first  super-rows  of  V  are  2x2  identity 
matrices.  Therefore,  in  the  first  super-row  of  V  there  are  three  (n)  2x2  (mxm)  identity 
matrices.  The  second  super-row  is  computed  using  equation  4.50  with  i  (the  super-row 
index}  equal  to  two  and  j  (the  super-column  index)  running  from  one  to  three  since 
there  are  three  super  columns  in  the  V  matrix. 


V    - 


1 

0 

1 

0 

1 

0 

0 

1 

0 

1 

0 

1 

xl 

0 

h 

0 

^3 

0 

0 

X, 

0 

*o 

0 

X 

(4.55) 


The  entries  of  matrix  G  are  given  by  equation  4.52.  Since  m=2  the  Lj  are  2x1  column 
vectors.  The  matrix  G  will  have  three  (n)  of  these  column  vectors  in  its  diagonal. 


G    = 


5n 

0 

0 

fei 

0 

0 

0 

<12 

0 

0 

^22 

0 

0 

0 

&13 

0 

0 

^23 

(4.56) 


If  one  multiplies  equation  4.55  (T)  by  equation  4.56  (£)  the  result  obtained  will 
be  equivalent  to  the  one  that  would  be  obtained  if  equation  4.47  is  used  for  the  same 
values  of  n,  m  and  €. 

In  order  to  find  the  initial  states  the  inversion  of  a  matrix  must  be  performed,  as 

T 

was  mentioned  at  the  beginning  of  this  Section.  The  matrix  R=QQ     is  the  one  that 

has  to  be  inverted  as  shown  in  equation  4.41.  Thus,  for  the  multiple-output  case  the 
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condition  number  of  R  is  the  value  that  gives  a  measure  of  how  observable  the  state 
variables  are.  Expanding  R  the  following  is  obtained: 

R    =    QQT  (4.57) 

=  (VGM)T(VGM) 
=  MTGTVTVG\l 

T 

In  general,  any  quadratic  form  of  the  type  x   Px(P  positive  definite)  is  bounded 

below  by  the  product  of  the  minimum  eigenvalue  of  P  and  the  length  of  x  and  bounded 
above  by  the  product  of  the  maximum  eigenvalue  of  P  and  the  length  of  x  [Ref.  13:  p. 
121].    For  the  particular  case  of  equation  4.57  this  result  can  be  written  as 


>-niinl|x||2  <  xTQQTx  <  Xmax||x||2  (4.58) 

where  Xmax  and  '^m^n  are  the  maximum  and  minimum  eigenvalues  of  the  square 
matrix  R=QQ   .  The  theorem  mentioned  in  [Ref.  13:  p.  121]  can  be  used  several  times 

in  equation  4.58  until  upper  and  lower  bounds  are  obtained  only  as  functions  of  the 

T  T  T 

maximum  and  minimum  eigenvalues  of  the  square  matrices  V  V,  G   G  and  MM. 

Equation  4.59  presents  what  happens  for  the  upper  bound 


xQQTx    =  xT\[TGT(VJV)G\[\  (4.59) 

-  >W^iixTmT(2t£)MxH 

^  Xmax(ZTZ)Xmax(GTG)||xTMTMx|| 
^Xmax(ZTZ)Xmax(GTG)Xmax(MTM)||x||2 

and  an  equivalent  result  can  be  obtained  for  the  lower  bound,  i.e., 

xQQTx>X^n(VTV)Xnmi(GTG)\mn(\lTM)\\x\\2  (4.60) 

Recalling  the  fact  that  in  general  the  norm  of  the  product  is  less  or  equal  than  the 
product  of  the  norms,  the  three  results  presented  in  equations  4.58.  4.59  and  4.60  imply 

that 

Xmax(QQT)^^max^T^max^Te)Xmax(MTM)  (4.61) 
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and 


In  [Ref.  16:  p.  38]  the  author  mentions  that  the  ratio  between  the  largest  and  smallest 
eigenvalues  is  a  measure  of  how  ill-conditioned  a  matrix  is.  With  this  in  mind  and 
looking  at  the  inequalities  4.61  and  4.62  it  can  be  said  that 

cond(QQT)  <  [cond(7Ty)][cond(GT£)][cond( MTM)]  (4.63) 

Observe  that  although  the  condition  number  of  QQ1  is  a  function  of  V*V,  G   G  and 

T 
MM.,  any  improvement  of  the  degree  of  observability  of  the  system  can  only  be 

T 

accomplished  by  changing  the  condition  number  of  G   G  because  this  is  the  only  part 

of  the  product   of  the  three  matrices   that   is   not  totally  dependent   on   the   plant. 

Assuming  that  the  user  is  only  allowed  to  alter  the  C  matrix  and  recalling  equation 

T 

4.43,  it  can  be  seen  that  only  by  minimizing  the  condition  number  of  G   G  will  the 

T  T      •  •  • 

condition  number  of  QQ     be  minimized.  The  matrix  G   G  is  a  nxn  diagonal  matrix 

whose  entries  are  actually  the  lengths  of  the  column  vectors  of  G  squared.  Let  the 

T 

entries  of  G   G  be  the  2:; 
-    -  ~]J 

§ll  =  I^il 


v —    2 

822  =  1  ^2 


i=l,2 m  (4.64) 


_  v"1  "    2 


p 


Thus,  the  condition  number  of  GTG  as  it  is  defined  in  [Ref.  16:  p.  38],  will  be  given  by 
the  ratio  of  the  largest  an 
is  given  by  equation  4.65: 


T  -    T 

the  ratio  of  the  largest  and  the  smallest  entries  in  GlG.    The  condition  number  of  G   G 


m 


I  Sit2 


=  (cond  GTG)  (4.65) 


m 

v  r.  2 

—  -is 
i=  l 
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where  (,■£  is  the  largest  element  in  the  G  G  matrix  while  i,-.  is  the  smallest  element  in 
the  same  matrix.  Comparing  equation  4.65  with  equation  4.34  and  recalling  the  fact 
that  the  numerator  and  denominator  of  equation  4.34  are  respectively  the  largest  and 
smallest  entries  in  the  matrix  G  (which  is  the  equivalent  to  the  matrix  G  for  the  multi- 
output  case),  it  can  be  seen  that  starting  from  two  different  definitions  of  condition 
number,  two  equivalent  equations  for  the  same  quantity  are  obtained.  The  entries  of 
the  matrix  G  are  also  summations  and  are  given  ,  in  general,  by 

n 

Sij  =  ^V^pj     i=1'2 m    J=l>2....,n  (4.66) 

p=l 

Substituting  equation  4.66  in  equation  4.65  for  the  values  of  Cp  and  C    yields 


m      n 


I  (Iw*' 


i=  l  p=  l  T 

=(condG1G)  (4.67) 

m      n 

I  dvV2 

1=  1  p=  1 


At  this  point,  if  one  attempts  to  use  the  same  method  used  in  the  single  output  case, 
namely  solving  equation  4.67  for  the  entry  of  C  that  is  related  to  the  state  whose 
relative  observability  is  to  be  improved,  it  can  be  seen  that  this  task  is  far  from  being 
an  easy  one  because  the  equations  that  have  to  be  solved  are  quadratics.  Therefore, 
one  of  the  possible  solutions  for  this  problem  is  the  following:  the  first  assumption  is 
that  the  user,  after  studying  the  possibilities  of  the  system  under  test,  knows  already 
which  entries  of  C  can  be  changed  and  how  much  can  they  be  changed.  Then  only  one 
of  the  entries  of  C  would  be  changed  in  each  run  over  the  allowed  interval.  This 

procedure  would  have  to  be  repeated  for  every  entry  in  the  C  matrix  that  can  be 

T 

changed.     Meanwhile,  using  equation  4.67  the  condition  number  of  Q  G  would  be 

T 
computed  as  well  as  the  condition  number  of  QQ    .  For  every  change  in  the  entry  ol  C 

under  test  the  relative  degree  of  observability  of  the  system  would  be  computed.  A  plot 

T  T 

of  the  condition  number  of  G  G  against  the  condition  number  of  QQ      for  every 
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possible  allowable  change  in  the  (."  matrix  would  be  an  extremely  valuable  output. 
With  this  information  and  the  knowledge  of  what  states  are  more  observable  under 
different  conditions,  a  choice  of  the  best  allowable  combination  of  the  entries  of  C 
could  then  be  made  such  that  the  condition  number  of  QQ1  is  minimum  or  a  specific 
state  variable  is  more  observable  than  the  others.  Notice  the  parallelism  that  exists 
between  the  derivations  involving  the  multi-output  case  and  the  generalized  graphical 
method  developed  in  Chapter  III  to  evaluate  the  relative  degree  of  observability  of  the 
state-variables.  In  both  cases  the  matrix  QQ  was  involved  as  the  center  of  the 
problem:  a  search  for  the  minimum  eigenvector  in  the  method  of  Chapter  III  and  the 
search  lor  the  minimum  condition  number  of  QQ  in  the  method  presented  in  this 
Chapter.  The  goals  of  both  methods  are  somehow  different  although  correlated 
because  the  second  uses  information  from  the  first  one. 

E.        CASE  WHERE  THERE  IS  A  CONTROL  INPUT 

In  the  beginning  of  Chapter  III  the  concept  of  relative  observability  was  defined 
as    a    measure    of  how   accurately    one    can    determine    the    initial    states    from   the 

measurements  of  the  output  at  time  t  =  ti,  t-> t  .    The  control  input  u(t)  was  set  to 

zero  because  the  effect  of  the  input  signal  can  always  be  determined  if  x(t(j  can  be 
found.  Now.  considering  the  existence  of  a  control  input  u(t)*0,  the  concept  of  relative 
observability  must  be  redefined  as  follows:    measuring  y(t)  and  u(t)  at  time  t   =    tj. 

t-> t     one  wants  to  be  able  to  evaluate  which  state  components  are  most  and  least 

accurately  determined  from  the  measurements  of  both  the  input  and  output.  It  will  be 
shown  in  this  section  that  the  fact  that  one  is  measuring  the  input  signal  as  well  as  the 
output  signal  doesn't  alter  in  any  way  what  has  been  shown  to  be  valid  for  the  case  of 
u(t)  =  0.  Consider,  again,  the  LTI  system  described  by  equations  2.1  and  2.2. 
Differentiating  equation  2.2  n-1  times  and  substituting  x  for  its  value  given  in  equation 
2.1.  the  following  is  obtained: 

y(t0)  =  Cx(t0)  (4.68) 

y(t0)=CAx(t0)  +  CBu(t0) 
y(t0)=CA2x(t0)+CABu(t0)  +  CBu(t0) 
y(t0)=CA3x(t0)  +  CA2Bu(t0)  +  CABu(t0)4-CBu(t0) 

fy(tn)  =  CAn-1x(tn)  +  CAn"2Bu(t0)  +  CAn-3Bu(t0)  +  ...  +  CBu(t0) 


:o 


Notice  that  the  arguments  of  the  input  and  its  derivatives  is  Iq  which  means  that  the 
derivatives  of  the  inputs  are  being  evaluated  at  time  t  =  tn.  Equation  4.68  can  he 
written  in  the  matrix  form 


yrft^-Q^ttf+nucrfto) 


(4.69) 


where: 

y^o)  is  a  nxl  column  vector  of  the  output  and  its  n-1  first  derivatives, 

x(tQ>  is  a  nxl  column  vector  of  the  state  variables  at  time  tQ 

n   is   a   nx(n-l)  lower  triangular  Toeplitz  matrix  that  represents  the  impulse 

response  of  the  system  and  is  given  by  equation  4.70 


n  = 


0 

0 

0 

CB 

0 

0 

CAB 

CB 

0 

CAn"2B 

CAn"3B 

c 

(4.70) 


and  u^(tQ)   is   a    lx(n-l)   column  vector   of  the  input   and  its   n-2   first   derivatives. 
Considering,  now,  the  single  output  case,  equation  4.69  is  solved  for  x(tQ) 


-/rJVL, 


TV1, 


x(t0)  =  (Q1)"1yd(t0)-(Q1)"1nud(t0) 


(4.71) 


Observe  that  the  projection  of  an  error  made  in  measurements  of  y(t)  and  u(t)  in 
computing  the  initial  states  is  determined,  again,  by  the  degree  of  ill-conditioning  of 
Q  which  is  the  only  matrix  that  needs  to  be  inverted  in  this  process.  The  matrix  Q 
causes  no  problem  to  because  it  can  be  determined  very  accurately  (assuming  that  the 
identification  of  the  plant  has  been  done).  Consider  the  example  used  in  Chapter  III 
to  explain  the  graphical  method  where 


Qt  = 


1.01   0.99 
0.99    1.01 


The  matrix  £1  was  equal  to  zero  because  the  control  input  was  set  to  zero.  Now  the 
control  input  is  different  from  zero  which  implies  that  £1  is  also  different  from  zero. 
Let  n  be  equal  to 
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0 
a 


where  alpha  is  any  real  number  different  from  zero.    Substituting  the  matrices  Q^  and 
n  in  equation  4.69  yields 

l.01x1(t0)+0.99x2(t0)  =  y1(t0)  (4.-2) 


0.99xj(t0)+  1.01x2(t0)  =  y2(t0)  -  aud(t0)  (4.73) 

Equations  4.72  and  4.73  can  be  plotted  in  the  Xj^q)  x2(t,0  plane  as  was  done  in 
Chapter  III  Section  E.  However,  just  by  inspection  one  can  see  that  the  slopes  of  the 
lines  corresponding  to  equations  4.72  and  4.73  are  the  same  as  the  slopes  of  the  lines 
corresponding  to  equations  3.20  and  3.21  although  a  new  member  exists  in  equation 
4.73.  This  new  member  appears  due  to  the  existence  of  a  control  input  and  is  only 
responsible  for  a  change  in  the  intercepts,  i.e..  for  a  change  in  the  solution  of  the 
system  of  equations  of  the  example  in  Chapter  III.  Since  the  slopes  are  equal  so  are 
the  angles  between  the  lines.  It  was  shown  in  Chapter  III  that  the  angles  between  these 

T 

lines  are  highly  correlated  with  the  degree  of  ill-conditioning  of  the  Q  matrix:  the 
smaller  the  angles  the  higher  the  degree  of  ill-conditioning  of  the  matrix  Q  .  Therefore. 
the  fact  that  a  control  input  exists  doesn't  bring  any  changes  in  the  basic  ideas 
presented  so  far  leading  to  the  evaluation  of  the  relative  degree  of  observability  of  the 
state  variables  of  a  given  system.    The  same  thing  happens  for  the  improvement  of  the 

degree   of  observability   of  a  weakly   observable   state.    In   the   case   of  u(t)=0   the 

T 
procedure  was  the  factorization  of  the  Q     matrix  as  the  product  of  three  full  rank 

matrices  such  that  only  one  was  dependent  on  the  C  matrix,  which  was  the  only  part 

of  the  system  that  the  user  is  allowed  to  change.  Then,  the  goal  was  to  find  a  new 

entry  in  the  C  matrix  such  that  the  condition  number  of  Q     was  minimum.  Once  a 

minimum  was  found  the  relative  observability  tests  were  performed  as  well  as  the 

T 

evaluation  of  the  new  condition  number  of  Q  .  In  this  case,  where  a  control  input 
exists,  the  same  procedure  is  followed  and.  furthermore,  the  results  obtained  will  be 
equal  to  the  ones  that  would  be  obtained  if  u(t)  =  0  because  the  observability  matrix  Q 
is  independent  of  u(t). 
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The  case  of  multi-output  with  control  input  is.  again,  very  similar  to  the  multi- 
output  case  with  u(t)  =  0.  Consider  equation  4.69  where  the  matrix  Q     is  no  longer  a 

T 

nxn  matrix.  In  order  to  compute  the  initial  states  the  direct  inversion  of  Q     cannot  be 

T 

made  because  Q  is  not  a  square  matrix  .  Premultiplying  equation  4.69  by  Q  equation 
4.74  is  obtained 

Qyd(t0)=QQTx(t())  +  Qnud(t0)  (4.74) 

and  solving  it  for  x(tQ) 

x(t0)  =  (QQT)'1yd(t0)-(QQ1V1nud(t0)  (4.75) 

T 

Observe  that  in  this  case  the  degree  of  ill-conditioning  of  QQ     is  the  one  that  has  to  be 

T 

studied  instead  of  the  degree  of   ill-conditioning  of   Q   .   For  the  evaluation  of  the 

relative  degree  of  observability  the  generalized  graphical  method  is  the  most  suitable 

t 

one  because  the  matrix  involved  in  the  derivation  of  the  algorithm  is  the  QQ     matrix. 

For  the  case  of  improvement  of  the  relative  observability  of  a  particular  state  variable 
the  problems  that  were  presented  in  the  Section  D  of  this  Chapter  are  also  valid  here 
because  the  square  matrix  QQ     is  independent  of  the  control  input. 

F.       SUMMARY 

In  this  Chapter  a  method  to  improve  the  relative  observability  of  a  weakly 
observable  state  was  developed  for  the  single  output  regulator  (u(t)  =  0)  case.  This 
method  was  based  on  the  concepts  of  spectral  factorization  of  a  matrix  with  distinct 

eigenvalues  discussed  in  [Ref.  10:  p.  609].    The  method  has  the  goal  of  minimizing  the 

T 

condition  number  of  the  Q     matrix.  An  algorithm  was  derived  and  a  Fortran  program 

(COBS)  presented  in  Appendix  B  was  written.  This  program  improves  the  relative 
degree  of  observability  of  the  least  observable  state.  Documentation  is  provided  to 
enable  the  user  to  improve  the  relative  degree  of  observability  of  a  particular  state- 
variable  other  than  the  least  observable  one.  This  method  was  then  extended  for  the 
multi-output  case  where  the  concept  of  the  pseudo-inverse  of  a  matrix  was  used.  The 
factorization  of  Q  into  the  product  of  three  matrices  where  only  one  was  not  entirely 
dependent  on  the  plant  was  performed  and  examples  were  presented  to  demonstrate 
that  this  general  case  factorization  contains  the  single  output  case.  The  method  to 
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improve  the  relative  degree  of  observability  of  a  particular  state-variable  was  then 
presented  and  shown  to  be  different  from  the  single  output  case  as  well  as  much  more 
elaborate.  Finally,  the  case  where  a  control  input  is  present  was  studied  both  for  the 
single  and  multi-output  cases,  leading  to  the  conclusion  that  the  fact  that  u(t)*0 
doesn't  change  in  any  way  what  has  been  said  for  the  cases  where  u(t)  =  0. 
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V.  ANALYSIS  OF  THE  RESULTS 


A.  INTRODUCTION 

In  the  last  two  Chapters  algorithms  have  been  developed  to  test  the  relative 
degree  of  observability  of  a  LTI  system  and  to  improve  the  degree  of  observability  of  a 
weakly  observable  state.  As  previously  mentioned,  these  algorithms  are  implemented  in 
two  programs  (TOBS  and  COBS)  whose  performance  will  be  discussed  in  this  Chapter. 
The  results  obtained  while  running  the  programs  will  be  evaluated  and  some 
observations  will  be  presented. 

B.  PERFORMANCE  OF  TOBS 

This  program  tests  and  ranks  the  relative  degree  of  observability  for  each  state 
variable  using  three  out  the  four  methods  presented  in  Chapter  III.  namely  the  upper- 
bound  error  method,  the  graphical  method  and  the  generalized  graphical  method.  In 
the  first  phase  of  testing  the  program,  the  performance  of  the  methods  was  investigated 
using  examples  already  tested  and  used  in  [Ref.  7].  The  Table  below  shows  the  values 
for  the  degrees  of  observability  obtained.  As  can  be  seen  all  the  three  methods  agree  on 
what  are  the  least  and  most  observable  states.  The  matrices  A  and  C  of  this  system 
are: 


A  = 


■12       3 

3        -4 


C  =      [0     1] 


TABLE  1 
TOBS  OUTPUT  FOR  ABLINS  SECOND  ORDER  EXAMPLE 


Upper-Bound 
Error 

Graphical 
Method 

Generalized 
Graphical 

x2 

1.0 

x2 

1.0 

x2 

0.68 

xl 

0.6 

xl 

0.3 

xl 

0.05 
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Table  1  shows  that  in  all  the  cases  x->  is  the  most  observable  state.  This  makes  sense 
because  the  matrix  C  shows  that  X2  is  the  state  variable  that  is  being  measured.  Once 
TOBS  was  tuned  and  running  as  desired,  several  state-space  representations  of  transfer 
functions  with  different  C  matrices  were  studied.  The  maximum  order  tested  was  order 
ten.  As  the  order  of  the  plant  starts  increasing  the  graphical  method  begins  disagreeing 
with  the  other  two  methods  in  the  ranking  of  the  least  observable  state.  However,  the 
results  of  the  upper-bound  and  generalized  graphical  methods  agree  with  one  another 
in  almost  all  cases.  This  behavior  can  be  understood  if  one  recalls  the  discussion  of  the 
graphical  method  presented  in  Chapter  III.  This  method  relies  on  finding  the  largest 
component  of  the  vector  relating  the  states  with  the  measurements.  As  the  order  of  the 
system  increases  the  surfaces  representing  the  relations  between  the  states  and  the 
measurements  start  being  planes  and  hyperplanes  where  the  graphical  reasoning  behind 
this  simple  procedure  is  no  longer  easy  to  apply.  However,  even  for  high  order  systems, 
the  most  observable  state  is  almost  always  the  same  in  all  three  methods.  The  big 
advantage  of  the  graphical  method  is  the  small  amount  of  computation  required, 
enabling  the  user  to  obtain  a  quick  determination  of  the  relative  observability  of  the 
system  in  low  order  or  reduced  order  models.  Furthermore,  due  to  its  relatively  easy 
graphical  interpretation,  it  can  be  quite  helpful  in  understanding  the  ideas  that  serve  as 
the  basis  for  this  type  of  method.  Although  agreeing  with  the  generalized  graphical 
method,  the  upper-bound  error  method  suffers  from  the  fact  that  it  needs  a  matrix 
inversion  to  obtain  the  required  results.  With  this  method,  the  evaluation  of  the 
relative  degree  o[  observability  of  the  state  variables  cannot  be  obtained  when  the 
system  is  not  completely  observable  because  the  inversion  of  the  observability  matrix  is 
needed.  In  any  of  the  other  methods  this  information  can  still  be  obtained  even  if  one 
or  more  states  are  not  observable.  The  generalized  method  seems  to  be  the  best  one 
because  although  it  presents  a  little  more  difficult  graphical  interpretation,  it  doesn't 
need  any  matrix  inversion  and,  more  importantly,  it  can  be  readily  used  to  compute  the 
relative  degree  of  observability  for  each  state  variable  in  the  multi-output  case  while 
the  other  two  methods,  as  they  are  developed,  need  some  changes  before  they  can  be 
used  in  evaluating  the  relative  degree  of  observability  for  each  state  variable  in  the 
multi-output  case. 

As  was  previously  mentioned,  several  runs  were  made  using  state-space 
representations  of  transfer  functions  from  order  two  to  order  ten  for  different 
combinations  of  the  entries  of  the  C  matrix.    Since  a  large  amount  of  computer  output 
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would  probably  not  be  of  much  help  in  understanding  the  relative  performance  of  the 
methods  described  in  Chapter  III.  a  summary  of  the  results  is  given  in  Table  2. 

TABLE  2 

SUMMARY  OF  THE  RESULTS  OBTAINED  USING  THREE  DIFFERENT 

METHODS 


Order  of 
the  Plant 

R12 

R13 

R23 

2 

100% 

100% 

100% 

3 

95% 

95% 

100% 

4 

93% 

93% 

99% 

5 

8S% 

87% 

97% 

6 

82% 

82% 

97% 

7 

75% 

73% 

9S% 

S 

71% 

71% 

96% 

9 

69% 

68% 

94% 

10 

65% 

65% 

94% 

In  Table  2  column  one  contains  the  order  of  the  plant.  Columns  two,  three  and 
four  contain  values  that  are  a  measure  of  the  correlation  between  the  graphical  method 
(number  one)  and  upper-bound  error  (number  two)  -R^"'  ^e  graphical  and  the 
generalized  graphical  method  (number  three)  -R13-  and  finally,  between  the  upper- 
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bound  error  method  and  the  generalized  graphical  method  -KjY-  The  values  in  Table  2 
were  computed  using  equation  5.1 

Rij    -  -^^  U- 1,2,3  (5.1) 

J  n  nr 

where  n  is  the  order  of  the  system  (number  of  states),  nr  is  the  number  of  runs 
performed  for  the  same  value  of  n  but  for  different  entries  of  the  C  matrix  and  nm  is 
the  number  of  matches  between  the  two  methods  insofar  as  the  relative  degree  of 
observability  is  concerned.  Notice  that  the  product  n  nr  represents  the  maximum 
possible  number  of  agreements.  For  instance,  suppose  that  the  order  of  the  system  is  8 
and  one  performs  10  runs.  The  maximum  possible  number  of  matches  is  SO.  However. 
if  the  two  methods  only  agree  in  60  (n  =60)  out  oi~  the  SO  possible,  the  agreement 
percentage  will  be  75%  and  this  is  the  value  that  appears  in  the  corresponding  columns 
of  Table  2.  Looking  at  this  Table  one  can  observe  what  was  mentioned  earlier  in  this 
Chapter,  namely  the  degradation  of  the  graphical  method  for  the  evaluation  of  the 
relative  degree  of  observability  as  the  order  of  the  system  increases.  Actually,  the 
disagreement  between  the  graphical  method  and  the  other  two  is  concentrated  in  the 
least  observable  states.  In  most  of  the  cases  the  most  observable  states  were  the  same 
within  one  level  of  ranking  in  all  three  methods.  As  can  be  seen  in  the  fourth  column 
of  Table  2.  the  agreement  percentage  is  slowly  decreasing  as  the  order  of  the  system 
increases.  This  occurs  due  the  way  the  values  are  computed.  In  fact,  the  values  for  the 
degrees  of  observability  for  the  most  observable  states  are,  in  most  cases,  of  the  same 
order  of  magnitude  for  high  order  systems.  A  slight  change  in  one  digit  can  produce  a 
different  ranking  in  the  state-variables,  which  is  interpreted  by  equation  5.1  as  a 
mismatch,  producing  a  decrease  in  the  agreement  percentage.  However,  as  the  order  of 
magnitude  of  the  degree  of  observability  is  the  same,  one  can  consider  that  the  two 
methods  still  agree  as  far  as  the  relative  degree  of  observability  is  concerned.  Of  course. 
this  decision  can  only  be  made  by  direct  analysis  of  the  computer  output.  To  illustrate 
this  characteristic  behavior  the  following  example  is  presented.  Consider  the  sixth  order 
svstem  where 
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A      = 


0 

1 

0 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

0 

1 

6 

-8 

-10 

-4 

-1 

-3 

c     = 


The  output  of  TOBS  for  this  system  is  given  in  Table  3. 

TABLE  3 
TOBS  OUTPUT  FOR  A  SIXTH  ORDER  EXAMPLE 


Upper-Bound 
Error 

Graphical 
Method 

Generalized 
Graphical 

x6 

1.0 

x5 

1.0 

x6 

0.999 

x5 

1.0 

x6 

1.0 

x5 

0.997 

x2 

0.66 

x2 

0.52 

x2 

0.915 

x4 

0.53 

x4 

0.08 

x4 

0.522 

xl 

0.51 

x3 

0.06 

x3 

0.421 

x3 

0.26 

xl 

0.05 

xl 

0.321 

As  can  be  seen  in  columns  four  and  six,  although  theoretically  the  upper-bound  error 
and  the  generalized  graphical  methods  do  not  agree  in  the  two  most  observable  states, 
they  are.  in  fact,  equally  observable  because  their  degree  of  observability  is  of  the  same 
order  of  magnitude  within  the  same  method.  This  just  shows  that  one  must  be  careful 
in  arriving  at  any  kind  of  conclusion  just  based  on  the  output  of  the  computer 
program. 
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The  third  and  probably  more  interesting  part  of  this  study  was  the  test  of  models 
of  real  plants.  The  first  model  studied  was  the  reduced  order  model  of  a  SL-7  hull  with 
three  degrees  of  freedom.  The  state  space  representation  was  derived  in  [Ref.  18:  pp. 
585.592]  and  is  shown  below  : 


v 

r 

,  T 


-2.64e-2 

-1.35e-4 
0 


-8.027  0 

■5.87e-2  0 

1  0 


r 


where: 


v  is  the  sway  rate 
r  is  the  yaw  rate 
¥  is  the  yaw 

Two  runs  of  this  model  were  performed  with  two  different  C  matrices.  The  results 
presented  in  Table  4.  show  the  output  of  TOBS  for  C  =  [  0    1    1  ] 

TABLE  4 
TOBS  OUTPUT  FOR  THE  SL-7  HULL  MODEL  -CASE  I- 


Upper-Bound 
Error 

Graphical 
Method 

Generalized 
Graphical 

x2 

1.0 

x2 

0.89 

x2 

0.99 

x3 

0.71 

x3 

0.39 

x3 

0.98 

xl 

0.02 

xl 

0.01 

xl 

0.01 

and  Table  5  shows  the  output  of  TOBS  for  C  =  [  0  0    1  ]. 

As  can  be  seen  in  Tables  4  and  5.  in  both  cases  the  three  tests  gave  v.  the  sway 
rate,  as  the  least  observable  state,  which  agrees  quite  strongly  with  one's  intuition 
about  the  problem.  Furthermore,  in  [Ref.  IS]  the  author  mentions  that  for  some  type 
of  studies  the  sway  rate  can  be  neglected  in  a  reduced  order  model  because  its  influence 
in  the  behaviour  of  the  third  order  model  is  not  verv  strons. 


65 


TABLE  5 
TOBS  OUTPUT  FOR  THE  SL-7  HULL  MODEL  -CASE  II- 


Upper-Bound 
Error 

Graphical 
Vie  mod 

Generalized 
Graphical 

x3 

1.0 

x3 

1.0 

x3 

1.0 

1.0 

x2 

1.0 

x2 

0.99 

xl 

0.01 

xl 

0.01 

xl 

0.01 

The  second  model  studied  in  this  phase  is  a  seventh  order  plant  that  is  presented 
as  example  two  of  [Ref.  19].  where  the  author  is  seeking  the  identification  of  the 
essential  states  to  reduce  the  order  of  the  model.  With  his  test  he  finds  that  states  x^  x^ 
x./:  and  \j  are  the  ones  that  contain  more  information  on  the  dominant  modes  of  the 
plant  and  are.  therefore,  the  ones  to  keep  in  the  reduced  order  model.  Performing  the 
three  tests  described  in  Chapter  III  on  his  model,  the  results  shown  in  Table  6  are 
obtained.  The  relative  observability  tests  performed  by  TOBS  show  that  the  most 
observable  states  are  x^,  x^.  x^  and  Xy  (with  special  strength  in  the  three  last  ones), 
which  agrees  with  the  results  obtained  by  the  authors  of  [Ref.  19].  An  interesting  point 
to  notice  in  this  example  is  that  there  is  a  difference  of  one  order  of  magnitude  in  the 
degrees  of  observability  of  x^,  Xi  and  x->  between  the  graphical  and  the  generalized 
graphical  method.  If  only  the  graphical  method  was  available  one  could  conclude  that 
x^.  x.,  x->  and  x->  were  all  weakly  observable.  However,  the  upper-bound  error  method 
gives  magnitudes  for  the  degrees  of  observability  of  the  state  variables  that  are  of  the 
same  order  of  magnitude  as  the  ones  given  by  the  generalized  graphical  method. 
Therefore,  in  this  example  only  x^,  and  perhaps  x?,  would  be  classified  as  weakly 
observable  states. 

C.       PERFORMANCE  OF  COBS 

The  algorithm  developed  in  Chapter  IV  for  the  single  output  case  lead  to  COBS. 
a  Fortran  program,  that  improves  the  relative  degree  of  observability  of  the  least 
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TABLE  6 
TOBS  OUTPUT  FOR  THE  MODEL  OF  EXAMPLE  2  OF  [REF.  19] 


Upper-Bound 
Error 

Graphical 
Method 

Generalized 
Graphical 

x6 

1.0 

x6 

1.0 

x6 

0.99 

x5 

1.0 

x5 

1.0 

x5 

0.99 

X- 

0.99 

x7 

0.79 

x7 

0.98 

x4 

0.69 

x4 

0.08 

x4 

0.92 

xl 

0.68 

xl 

0.05 

xl 

0.81 

x2 

0.11 

x2 

0.04 

x2 

0.80 

x3 

0.01 

x3 

0.01 

x3 

0.02 

observable  state.  Several  runs  were  made  of  plants  of  order  two  up  to  order  seven.  For 
a  few  cases  of  the  seventh  order  plants  the  program  took  a  long  time  to  perform  the 

T  .... 

minimization  of  the  condition  number  of  Q  .  By  relaxing  the  minimization  stopping 
criterion  the  computing  time  decreased  a  bit  but  not  dramatically.  However,  the 
results  that  were  obtained  were  very  satisfactory  because  the  least  observable  state 
became  the  most  or  the  second  most  observable  state  in  all  the  cases  with  the 
corresponding  minimization  of  the  condition  number  of  Q    . 

An  interesting  aspect  to  observe  is  what  happens  before  and  after  the 
minimization  from  the  graphical  point  of  view.  Consider,  as  an  example,  the  second 
order  system  presented  in  [Ref.  7]  where  initially 


A  = 


12 


C  =      [0     1] 
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and  the  least  observable  state  is  Xj  as  shown  in  Table  1.  Before  the  minimization  the 

T 

condition  number  of  Q     is  11.667  where 


J 


1 

-4 


After  running  COBS  the  new  C  matrix  is  C  =  [  1.17  1  ].  Notice  that  the  entry  c,, 
was  the  one  changed  during  minimization  because  this  entry  corresponds  to  the  least 
observable  state  whose  relative  observability  is  to  be  improved.  With  this  new  C  matrix 
the  transpose  of  the  observability  matrix  is 


Qt  = 


1.17 
11.00 


1.0 
-0.5 


whose  condition  number  is  now  equal  to  2.S9.  The  new  ranking  of  the  state  variables 
based  on  the  relative  degree  of  observability  is  shown  in  Table  7.  Notice  again  the 
difference  of  one  order  of  magnitude  between  the  graphical  and  generalized  graphical 
methods.  The  problem  is  now  the  inverse  of  what  it  was,  i.e.,  now  the  graphical 
method  agrees  with  the  upper-bound  error  method  in  the  order  of  magnitude  of  the 
degrees  of  observability. 

TABLE  7 
NEW  RANKING  OF  THE  STATE  VARIABLES  FOR  ABLINS  EXAMPLE 


Upper-Bound 

Error 

Graphical 
Method 

Generalized 
Graphical 

xl 

0.91 

xl 

0.85 

xl 

0.84 

x2 

0.65 

x2 

0.53 

x2 

0.01 

The  goal  was  reached,  i.e.,  the  state  that  once  was  the  least  observable  one    is 

T 

now  the  most  observable  and  the  condition  number  of  Q     decreased. 

In  order  to  obtain  a  graphical  interpretation  of  what  occurs  in  the  system  by 
changing  the  C  matrix,  the  procedure  used  in  Chapter  III  Section  E  can  be  used  again. 
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Substituting  the  Q     matrix  before  the  minimization  in  equation  3.19  the  equations  of 

two  lines  will  be  obtained.  These  two  lines  can  be  plotted  in  the  Xj(t()).  x->(t,.)  plane. 
assuming  that  Vj  and  y^  are  given  measured  values.  For  plotting  purposes  y,  and  y-> 
are  considered  proportional  to  the  two  horizontal  lines  presented  at  the  top  of  Figures 

5.1  and  5.2.  Looking  at  Figure  5.1  it  can  be  seen  that  the  uncertainty  region  for  x->(t()) 
is  equal  to  the  uncertainty  region  produced  by  the  errors  in  the  measurements.  If  one 
relates  these  results  with  those  of  Table  1  one  can  see  that  they  agree  because  Xt  is  the 
most  observable  state.  Furthermore,  two  methods  give  for  this  state  variable  a  value  of 
1  for  the  degree  of  observability  which  means  that  an  error  in  the  measurements  is 
projected  with  the  same  order  of  magnitude  when  the  computation  of  the  initial  states 
takes  place.  The  state-variable  Xj  is  the  least  observable  because  its  uncertainty  region 
is  larger  than  the  one  produced  by  an  error  in  the  measurements.  Substituting  Q*  after 
minimization  in  equation  3.19  another  different  set  of  equations  relating  the  states  and 
the  measurements  of  the  output  are  obtained.    These  two  lines  are  presented  in  Figure 

5.2  in  the  usual  Xj(t(j).  x->(t(y)  plane.  The  distances  that  correspond  to  yj  and  y->  are 
now  the  double  of  what  they  were  in  Figure  5.1  but  their  ratio  is  the  same.  This 
rescaling  was  done  to  obtain  a  cleaner  plot  of  the  two  lines.  Notice  that  now  the 
uncertainty  region  for  *ia(t())  is  of  the  same  order  of  magnitude  as  the  one  produced 

by  the  measurements  of  the  output  while  the  uncertainty  region  for  x->a(tQ)  increased. 

T 

The  condition  number  of  Q     decreased  because  the  smaller  angle  between  the  lines  in 

Figure  5.1  is  about  33°  and  in  Figure  5.2  is  about  42°. 

In  the  runs  of  COBS  with  the  different  plants  in  about  40%  of  the  cases  the 
initially  least  observable  state  became  the  second  most  observable,  but  with  a  value  for 
the  degree  of  observability  of  the  same  order  of  magnitude  as  the  most  observable 
state.  Unfortunately,  in  the  real  world,  the  user  may  not  be  able  to  change  the  C 
matrix  according  to  the  least  observable  state.  Generally  there  are  constraints  both  for 
the  entries  of  the  C  matrix  that  one  is  allowed  to  change  as  well  as  for  the  values  that 
the  entries  can  have.  The  program  COBS  allows  a  change  in  the  entry  of  the  C  matrix 
chosen  by  the  user  although  the  default  is  the  the  change  in  the  entry  of  C  associated 
with  the  least  observable  state.  As  previously  mentioned,  documentation  in  the  listing 
of  the  program  is  provided  to  guide  the  user  to  set  the  program  to  change  any  other 
entry  in  C.  Despite  the  somewhat  unrealistic  philosophy  upon  which  this  algorithm  was 
based,  it  provides  a  way  of  checking  how  the  system  and  the  states  behave  with  respect 
to  one  another,  which  can  be  helpful  to  earn'  out  further  investigations. 
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Figure  5.1     Graphical  interpretation  of  Ablins  example 
before  improvement. 
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Figure  5.2     Graphical  interpretation  of  Ablin's  example 
alter  improvement. 
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D.       SUMMARY 

In  this  Chapter  the  performance  of  the  methods  used  to  evaluate  the  relative 
degree  of  observability  of  a  set  o[  state  variables  was  analyzed  and  some  important 
results  were  obtained  from  the  use  of  the  program  that  performs  these  tests  (TOBS). 
The  analysis  of  the  performance  of  the  algorithm  in  improving  the  relative  degree  of 
observability  of  a  weakly  observable  state  was  done  and  the  graphical  interpretation  of 
a  second  order  plant  before  and  after  the  improvement  was  studied.  Finally  the 
inherent  advantages  and  disadvantages  of  this  algorithm  were  pointed  out. 
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VI.  CONCLUSION 

This  thesis  further  develops  the  issue  raised  by  Brown  [Ref.  2]  concerning  the 
question  of  "how  observable  a  linear  system  is".  This  issue,  which  leads  to  the  concept 
of  relative  observability,  is  addressed  in  four  different  ways  in  the  thesis.  Two  of  the 
methods  used  to  test  the  relative  observability  of  an  LTI  system  are  based  on  graphical 
interpretations  of  the  relations  that  exist  between  the  measurements  of  the  output  and 
the  states.  In  the  first  graphical  method  the  information  imbedded  in  the  angles 
between  the  vectors  whose  components  are  the  row  entries  of  the  observability  matrix 
and  the  coordinate  axes  (xj(tQ),  X2(tQ),...,  xn(t0))  is  used  to  determine  the  relative 
observability  of  the  state  variables.  The  smaller  the  angle  between  any  vector  and  an 
axis  of  the  state-space,  the  smaller  the  uncertainty  region  and  the  higher  the  degree  of 
observability  for  that  state-variable.  Since  the  angles  between  the  vectors  and  the  axes 
are  directly  related  to  the  magnitudes  of  their  components  along  these  axes,  an  easy 
algorithm  to  check  the  relative  degree  of  observability  of  a  state-variable  was  derived, 
constituting  one  of  the  three  methods  contained  in  the  computer  program  TOBS. 
During  the  study  of  this  method  it  was  observed  that  a  weakly  observable  state  was 
more  sensitive  to  changes  in  the  magnitudes  of  the  errors  in  the  measurements  than  a 
strongly  observable  state.  The  second  graphical  method  is  called  the  generalized 
graphical  method  and  is  developed  from  the  idea  of  searching  for  the  direction  of 
maximum  error  yielding  a  minimum  for  the  observability  function.  This  direction  turns 

T 

out  to  be  given  by  the  eigenvector  associated  with  the  smallest  eigenvalue  of  the  QQ 

matrix.    Therefore,   the  smallest  component  of  the  eigenvector  associated  with  the 

T 

smallest  eigenvalue  of  QQ     corresponds  to  the  most  observable  state  while  the  largest 

component  corresponds  to  the  least  observable  state. 

Two  other  analytical  methods  were  presented  and  explained:  the  upper-bound 
error  method  and  the  standard  deviation  error  method.    Although  these  two  methods 

seem  to  be  conceptually  different  from  the  two  graphical  methods,  a  common  thread 

T 

was  the  fact  that  the  degree  of  ill-conditioning  of  the  Q     matrix  plays  an  important 

role  in  the  theoretical  development  of  both  graphical  and  analytical  methods  because  a 
highly  ill-conditioned  matrix  has  a  twofold  meaning.  On  one  hand  it  means  that  the 
rows  of  that  particular  matrix  have  a  high  degree  of  dependency  (nearly  parallel)  and 
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on  the  other  hand  it  means  that  the  determinant  of  that  matrix  is  a  very'  small  number 
(which  causes  large  entries  in  the  inverse  of  the  matrix).  The  high  degree  of  dependency 
among  the  rows  of  Q1  is  one  oi"  the  basic  ideas  behind  the  development  of  the 
graphical  method  while  the  small  number  for  the  determinant  is  one  of  the  key  ideas 
behind  the  understanding  of  both  the  upper-bound  error  and  the  standard  deviation 
error  methods.  The  graphical,  upper-bound  error  and  generalized  graphical  methods 
were  throughly  tested  in  TOBS  using  different  models.  All  three  methods  agreed  in 
almost  all  cases  in  ranking  the  relative  observability  of  the  states.  However,  the 
graphical  method  degraded  as  the  order  of  the  system  increased.  Despite  this  feature, 
the  graphical  method  has  a  big  advantage:  it  requires  a  very  small  amount  of 
computation,  enabling  the  user  to  obtain  a  quick  determination  of  the  relative 
observability  in  low  order  and  reduced  order  system  models.  Furthermore,  due  to  its 
relatively  easy  graphical  interpretation  it  can  be  very  helpful  in  understanding  both  the 
concept  of  an  ill-conditioned  matrix  and  the  role  that  this  characteristic  of  a  specific 
matrix  has  in  computing  the  relative  observability  of  a  system.  The  upper  bound  error 
method  is  reliable  but  has  two  disadvantages:  the  requirement  of  matrix  inversion  for 
ever}'  evaluation  of  the  degree  of  observability  and  the  impossibility  of  using  the 
method  when  the  system  is  not  completly  observable.  The  generalized  graphical 
method  is  the  approach  whose  performance  is  more  satisfactory  because  it  can  be  used 
even  if  the  system  is  not  completly  observable,  doesn't  degrade  when  the  order  of  the 
system  increases  and  doesn't  need  any  matrix  inversion  for  the  required  results  to  be 
obtained.  Furthermore,  it  can  be  readily  used  to  compute  the  relative  degree  of 
observability  of  a  multi-output  system  while  the  graphical  and  upper-bound  error 
methods  require  some  changes  before  being  ready  to  compute  the  relative  degree  of 
observability  for  each  state-variable. 

In  the  second  part  of  this  thesis  an  algorithm  to  improve  the  relative  degree  of 
observability  for  the  single  output  regulator  case  was  developed  and  tested  in  the 

computer  program  COBS.   This  algorithm  was  based   on  the   concepts   of  spectral 

T 
factorization  of  a  matrix,  allowing  the  decomposition  of  the  Q     matrix  as  the  product 

of  three  matrices,  two  of  them  being  fully  dependent  on  the  the  plant  and  the  other 

depending  on  both  the  plant  and  the  C  matrix.  By  altering  the  entry  in  the  C  matrix 

associated  with  the  least  observable  state,  the  minimization  of  the  condition  number  of 

Q     was  achieved  and  the  least  observable  state  became  the  most  or  the  second  most 

observable  state  in  all  the  cases  that  were  tested.  The  single-output  regulator  case  was 
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then  generalized,  first  to  a  multi-output  regulator  problem  and  finally  to  a  multi-output 
servo  problem.  The  use  of  the  concept  of  the  matrix  pseudo-inverse  was  needed  for  the 
theoretical  development  o['  the  multi-output  problem  and  the  spectral  factorization  of 
Q  was  found  to  be  possible  although  hard  to  accomplish.  Two  illustrative  examples 
were  presented.  An  algorithm  to  improve  the  relative  degree  of  observability  of  a 
specific  state  was  also  proposed,  assuming  the  previous  knowledge  of  what  entries  in 
the  C  matrix  could  be  changed  and  how  much  could  they  be  altered.  Finally  it  was 
shown  that  the  servo  problem  added  nothing  new  to  what  had  been  accomplished  for 

T 

both  the  single  output  and  the  multi-output  case  because  the  key  matrices  (Q     for  the 

T 

single  output  and  QQ1    for  the  multi-output)  are  always  independent  of  the  control 

input. 

In  conclusion,  amongst  the  methods  developed  to  test  the  relative  degree  of 
observability  of  the  system,  it  seems  that  the  generalized  graphical  method  is  the  most 
reliable  one  although  the  graphical  method  can  give  quick  and  easily  computed 
information  about  the  most  observable  states.  Since  the  upper-bound  error  method 
needs  a  matrix  inversion,  it  may  become  computationally  involved  for  high  order 
systems  and  therefore  its  use  should  be  considered  depending  on  the  order  of  the 
system  to  be  studied. 

By  using  the  algorithms  developed  in  this  thesis,  the  designer  can  obtain  much 
insight  into  a  complex  system  with  a  relatively  small  amount  of  computations.  Exactly 
how  these  ideas  may  be  used  will  depend  on  the  designer's  knowledge  of  the 
capabilities  of  the  system  as  well  as  the  specifications  of  the  problem  under  study. 

There  are  two  recommendations  for  further  study.  The  first  one  is  the 
implementation  and  test  of  the  algorithm  proposed  for  the  improvement  of  the  degree 
of  observability  of  the  system  for  the  multi-output  case.  The  second  one  is  the  study  of 
the  issues  addressed  in  this  thesis  for  linear  time-varvins  svstems. 
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APPENDIX  A 
TESTING  THE  DEGREE  OF  OBSERVABILITY  OF  AN  LTI  SYSTEM 


PROGRAM  TOBS 
C 

c 

£*******************  PURPOSE  **************************** 

C  THIS  PROGRAM  EVALUATES  THE  DEGREE  OF  OBSERVABILITY  OF  * 

C  A  SYSTEM,  WHOSE  MATRIX  A  AND  C  ARE  IN  A  DATA  FILE,  BY  * 

C  THREE  DIFFERENT  METHODS:  A  GRAPHICAL  METHOD,   UPPER  * 

C  BOUND  ERROR  METHOD  C  AND  THE  GENERALIZED  GRAPHICAL  * 

C  METHOD.  THE  MAX  C  ORDER  FOR  THE  PLANT  IS  ORDER  TEN.  * 

C  THE  OUTPUT  IS  GIVEN  C  IN  A  DATA  FILE  WHERE  THE  STATES  * 

C  ARE  PRESENTED  IN  DECREASING  ORDER  OF  RELATIVE  OBSERV.  * 

C  INPUT  DATA  IS  FILE  FT04F001  -  LOGICAL  RECORD  IS  80  * 

C  OUTPUT  DATA  IS  FILE  FT08F001  -  LOGICAL  RECORD  IS  133  * 

c 
c 

c*****  DECLARATION  OF  VARIABLES  ***** 
C 

IMPLICIT  REALM  (A-H,0-X) 

DIMENSION  A(50,50) ,QT(50,50) , STATE(50) ,TEST(50) ,DGOBS(50) 

DIMENSION  DGTEST(50),C(1,50) ,QTN(50,50) ,QTNI(50,50) ,FLAG(50) 

DIMENSION  WKAREA(130) ,QN(50,50) ,QQN(50,50) 

COMPLEX  W(50) ,Z(50,50) 
C 
C 

c*****  READ  IN  ORDER  OF  SYSTEM  N,  AND  MATRICES  A  AND  C  ****** 
C 

READ(4,*)N 
C         WRITE(6. 600)1 

READ(4,*)(C(1,I),I=1,N) 

DO  20  1=1, N 
C  WRITE(6.700)I,J 

READ(4,*)(A(I,J),J=1,N) 
20    CONTINUE 
C 

c 

c*****  OUTPUT  MATRICES  A  AND  C  ***** 
C 

WRITE(8,450) 

WRITE (8, 350 )(C( 1,1), 1=1, N) 

WRITE(8,550) 

DO  25  1=1, N 


WRITE(8,350)(A(I,J),J=1,N) 


25 
C 

CONTINUE 

C 

C*****  COMPUTE  Q  TRANSPOSE  ****** 

n 

30 
C 

40 
50 
60 
C 

DO  30  J=1,N 

QT(1,J)=C(1,J) 
CONTINUE 

DO  60  1=2, N 

DO  50  JL=1,N 

QT(I,JL)=0. 
DO  40  JC=1,N 

OT(I,JL)=QT(I,JL)+QT(I-l,JC)*A(JC,JL) 
CONTINUE 
CONTINUE 
CONTINUE 
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c 

c*****  NORMALIZE  Q  TRANSPOSE  BY  ROW  ***** 
C 

DO  80  1  =  1, N 
SUM=0. 
DO  70  J=1,N 

SUM=SUM+QT ( I , J ) *QT ( I , J ) 
70        CONTINUE 

IF(SUM.EQ.O.)THEN 

QTN(I,JC)=QT(I,JC; 
ON(JC,I)=QTN(I,JC 

ELSE 

DO  75  JC=1,N 

QTN(I,JC)=QT(I,JC)/SQRT(SUM) 
QN(JC,I)=QTN(l,JC) 
75  CONTINUE 

END  IF 
80    CONTINUE 
C 

c 

C*****  FIND  BIGGEST  ENTRY  IN  EACH  COLUMN  ***** 
C 

DO  100  J=1,N 
BIG=0. 
DO  90  1  =  1, N 

IF(ABS(QTN(I, J)) .GT.BIG)THEN 

BIG=ABS(QTN(I,J)) 
END  IF 
90       CONTINUE 

FLAG(J)=1. 
IF(BIG.EQ.O. )THEN 
WRITE(6,800)J 
WRITE(8,800)J 
FLAG(J)=0. 
END  IF 

STATE(J)=BIG 
TEST(J)=STATE(J) 
100   CONTINUE 
C 
C 

C*****  COMPUTE  THE  RELATIVE  DEGREE  OF  OBSERVABILITY  ****** 
c*****  USING  THE  UPPER-BOUND  ERROR  METHOD  ****** 
C 

NSTOP=l. 

DO  105  1=1, N 

NSTOP=NSTOP*FLAG(I) 
105   CONTINUE 

IF(NSTOP.EQ.0.)THEN 

WRITE (6, 850) 

WRITE(8,850) 

END  IF 

IF(NSTOP.EQ.l.)THEN 

CALL  LINV2F(QTN,N,50,QTNI,0,WKAREA,IER) 
DO  120  1=1, N 
DEN=0. 
DO  110  J=1,N 

DEN=DEN+ABS ( QTNI ( I , J ) ) 
110  CONTINUE 

DGOBS(I)=l./DEN 
DGTEST(I)=DGOBS(I) 
120       CONTINUE 
END  IF 
C 
C 

c*****  MULTIPLYING  QN  BY  QTN  ***** 
C 

DO  126  1=1, N 

D0124  JL=1,N 

QQN(I,JL)=0. 
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DO  122  JC=1,N 

QQN ( I , JL ) =QQN ( I , JL ) +QN ( I , JC ) *QTN ( JC , JL ) 
122  CONTINUE 

124       CONTINUE 
126   CONTINUE 


THE  SMALLEST  EIGENVALUE 
EIGENVECTOR 


AND 


c*****  FINDING 

Qkkkkkkkkkkk 

c 

CALL  E I GRF ( QQN  ,N,50,1,W,Z,50, WKARE A , I ER ) 

SNALL=1000. 

DO  200  1  =  1, N 

IF(REAL(W(I)) .LT. SMALL) THEN 
SMALL=REAL(W(I)) 
ISMALL=I 
END  IF 
200   CONTINUE 


CORRESPONDING  ***** 

kkkkkkkkkkk 


c 
c 

Qkkkkk 

PRINT  IN  DECREASING  ORDER  THE  DEGREE  OF  OBSE 

Qkkkkk 
C 

FOR  EACH  STATE  VARIABLE 

WRITE(8,880) 

DO  140  K=1,N 

BPRIN1=0. 

BPRIN2=0. 

BPRIN3=0. 

DO  130  1=1, N 

IF  (TEST(I) .GE.BPRIN1)THEN 

BPRIN1=TEST(I) 

IBIG1=I 

END  IF 

IF  (DGTEST(I) .GT.BPRIN2.AND.NSTOP.EQ 

BFRIN2=DGTEST(I) 

IBIG2=I 

END  IF 

QTEST= ( 1 . - ABS ( REAL ( Z ( I , I SMALL ) ) ) ) 
IF(QTEST.GE.BPRIN3)THEN 

BPRIN3=QTEST 

IBIG3=I 

END  IF 

130 
C 

CONTINUE 

c 

Qkkkkk 
C 

kkkkkkkkkkk    OUTPUT  kkkkkkkkkkkkkkkkkk 

kkkkk 
kkkkk 


l.)THEN 


IF(NSTOP.EQ.0.)THEN 

WRITE(8,900)IBIG1,BPRIN1,IBIG3,BPRIN3 
END  IF 
IF(NSTOP.EQ.l.)THEN 

WRITE(8,910)IBIG1,BPRIN1,IBIG2,BPRIN2,IBIG3,BPRIN3 
END  IF 


140 


TEST(IBIG1)=-1. 
DGTEST(IBIG2)=-1. 
Z(IBIG3,ISMALL)=CMPLX(' 
CONTINUE 


0.) 


kkkkk    FORMATS  ***** 


350 
370 
450 
550 
560 
570 
800 
850 


FORMAT i 
FORMAT' 
FORMAT^ 
FORMAT i 
FORMAT I 
FORMAT i 
FORMAT^ 
FORMAT i 


\F13 
',F13 

ABLIN 


,7(2X,F13.7 
,7(4X,F13.7 
THE  ENTRIES 
THE  ENTRIES 
Q  TRANSPOSE 
Q  TRANSPOSE 


OF 
OF 
IS 


ARE 
ARE 


MATRIX 
MATRIX 

:') 

NORMALIZED  IS  : ' ) 
STATE  X1, 12, IX, 'IS  NOT  OBSERVABLE1,//) 
,'S  CRITERION  CAN  NOT  BE  EVALUATED,  ',/, 


BECAUSE 


7S 


&Q  TRANSPOSE  IS  SINGULAR',/) 
830   FORMAT ('  ',/,'  RELATIVE  DEGREE  OF  OBS .  FOR  GRAPHICAL  ABLIM1,1  AND 

&  GENERALIZED  METHODS1) 
900   FORMAT ( '  ' , 2X, ' X ' , 12 , 2X, F9 . 6 , 2X, 12 ( ' . ' ) , 2X, ' X ' , 12 , F9 . 6) 
910   FORMAT ( '  ' , 2X, ' X1 , 12 , 2X, F9 . 6 , 2X, ' X ' , 12 , 2X, F9 . 6 , 2X, ' X ' , 12 , 2X, F9 . 6 ) 
END 
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APPENDIX  B 
IMPROVING  THE  DEC  OF  OBS.  OF  A  WEAKLY  OBSERVABLE  SYSTEM 

PROGRAM  COBS 
C 

c 

Qkkkkkkkkkkkkkkkkkkkkkk     PURPOSE  ********************************** 

C  THIS  PROGRAM  COMPUTES  THE  RELATIVE  DEGREE  OF  OBSERVABILITY  * 

C  OF  THE  STATE  VARIABLES  OF  A  GIVEN  SYSTEM  BY  THREE  DIFFERENT  * 

C  METHODS  AND  THEN  COMPUTES  THE  NEW  ENTRY  IN  THE  C  MATRIX  IN  * 

C  ORDER  FOR  THE  DEGREE  OF  ILL-CONDITIONING  OF  THE  OBSERVABILITY  * 

C  MATRIX  TO  BE  MINIMUM.  IN  THIS  PROGRAM  THE  ENTRY  OF  C  TO  BE  * 
C  CHANGED  IS  THE  ONE  CORRESPONDING  TO  THE  LEAST  OBSERVABLE  STATE.* 

C  IF  THE  USER  WANTS  OR  NEEDS  TO  CHANGE  OTHER  ENTRY  OF  C  RATHER  * 

C  THAN  THE  ONE  ASSOCIATED  WITH  THE  LEAST  OBSERVABLE  STATE,  THE  * 

C  INTEGER  VARIABLE  LR  NEEDS  TO  BE  SET  FOR  THE  ORDER  OF  THE  STATE  * 

C  TO  BE  CHANGED.  FOR  INSTANCE,  IF  X4  NEEDS  TO  BE  CHANGED  INSTEAD  * 

C  OF  THE  ENTRY  ASSOCIATED  WITH  THE  LEAST  OBSERVABLE  STATE,  THE  * 

C  USER  SHOULD  PUT  LR=4  * 

Q     kkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkk 

c 
c 

c*****  DECLARATION  OF  VARIABLES  ***** 
C 

IMPLICIT  REAL*4  (A-H,0-X) 

DIMENSION  A(10,10) ,QT(10,10) , STATE ( 10) .TEST (10) ,DGOBS(10) 

DIMENSION  DGTEST(IO) ,C(1,10) ,QTN(10,10) ,QTNI (10 , 10) ,FLAG(10) 

DIMENSION  WK0(130).WK1(130) ,WK2(130) 

DIMENSION  QN(10,10) ,QON(10.10) , P (10 , 10) ,M( 10 , 10) 

DIMENSION  V(10,10),G{10,10) ,DA(10 , 10) ,DAR( 10 , 10) ,QTI (10 , 10) 

DIMENSION  AL(1,100) ,DARI(10,10) 

COMPLEX  W(10) ,Z(10,10) 
C 
C 

c*****  READ  IN  ORDER  OF  SYSTEM  N,  AND  MATRICES  A  AND  C  ****** 
C 

NT=0 

DEL=0.001 

READ(4,*)N 

READ(4,*)(C(1,I) ,1=1, N) 

DO  20  1=1, N 

READ(4,*)(A(I,J),J=1,N) 
20    CONTINUE 
C 
C 

c*****  OUTPUT  MATRICES  A  AND  C  ***** 
C 

WRITE(8,450) 

WRITE ( 8, 350 )(C (1,1), 1=1, N) 

WRITE(8,550) 

DO  25  1  =  1, N 

WRITE (8, 350) (A (I, J) ,J=1,N) 
25    CONTINUE 
C 
C 

Qkkkk-k   COMPUTE  Q  TRANSPOSE  ****** 
C 

CALL  MULT1(C,A,N,QT) 
C 
C 

c*****  NORMALIZE  Q  TRANSPOSE  BY  ROW  ***** 
C 

WRITE(6,700) 
700   FORMAT ('  ','  BEGINING  COMPUTING  THE  DEGREES  OF  OBSERVABILITY') 

80 


26    CONTINUE 

DO  SO  1=1, N 
SUH=0 . 
DO  70  J=1,N 

SUM=SUM+QT ( I , J ) *QT ( I , J ) 
70        CONTINUE 

IF(SUM.EQ.O.)THEN 

QTN(I,JC)=QT(I,JC) 
QN(JC,I)=QTN(I,JC) 
ELSE 

DO  75  JC=1,N 

QTN(I,JC)=QT(I,JC)/SQRT(SUM) 

QN(JC,I)=QTN(l,JC 
75  CONTINUE 

END  IF 
80    CONTINUE 
C 
C 

c*****  FIND  LARGEST  ENTRY  IN  EACH  COLUMN  ***** 
C 

DO  100  J=1,N 
BIG=0. 
DO  90  1=1 ,N 

IF(ABS(QTN(I, J)) .GT.BIG)THEN 

BIG=ABS(QTN(I,J)) 
END  IF 
90       CONTINUE 

FLAG(J)=1. 
IF(BIG.EQ.O.)THEN 
WRITE(6,800)J 
WRITE (8, 300) J 
FLAG(J)=0. 
END  IF 

STATE(J)=BIG 
TEST(J)=STATE(J) 
100   CONTINUE 
C 

c 

C*****  COMPUTE  THE  RELATIVE  DEGREE  OF  OBSERVABILITY  ****** 
c*****  USING  THE  UPPER-BOUND  ERROR  METHOD  ****** 
C 

NSTOP=l  . 

DO  105  1=1, N 

NSTOP=NSTOP*FLAG(I) 
105   CONTINUE 

IF(NSTOP.EO.O. )THEN 

WRITE(6,850) 

WRITE(3,850) 

END  IF 

IF(NSTOP.EQ.l.)THEN 

CALL  LINV2F(QTN,N,10,QTNI,2,WK0,IER) 
DO  120  1  =  1, N 
DEN=0 . 
DO  110  J=1,N 

DEN=DEN+ABS ( QTNI ( I , J ) ) 
110  CONTINUE 

DGOBS(I)=l./DEN 
DGTEST(I)=DGOBS(I) 
120      CONTINUE 
END  IF 
C 
C 

C*****  MULTIPLYING  QN  BY  QTN  ***** 
C 

DO  126  1  =  1, N 

DO  124  JL=1,N 

QQN(I,JL)=0. 
DO  122  JC=1,N 

QQN(I,JL)=QQN(I,JL)+ON(I,JC)*QTN(JC,JL) 

SI 


122  CONTINUE 

124       CONTINUE 

126   CONTINUE 
C 
C 
c*****  FINDING  THE  SMALLEST  EIGENVALUE  AND  CORRESPONDING  ***** 

Qkkkkkkkkkkk  EIGENVECTOR  kkkkkkkkkkk 

C 

CALL  EIGRF(QQN,N,10,1,W,Z,10,WK0,IER) 

SNALL=1000. 

DO  200  1  =  1, N 

IF(REAL(W(I)) .LT. SMALL) THEN 
SMALL=REAL(W(I)) 
ISMALL=I 
END  IF 
200   CONTINUE 
C 
C 

C*****  PRINT  IN  DECREASING  ORDER  THE  DEGREE  OF  OBSERVABILITY  ***** 
C*****  FOR  EACH  STATE  VARIABLE  kkkkk 

C 

WRITE (8, 880) 
DO  140  K=1,N 
BPRIN1=0. 
BPRIN2=0. 
BPRIN3=0. 
DO  130  1  =  1, N 

IF  (TEST(I) .GE.BPRIN1)THEN 
BPRIN1=TEST(I) 
IBIG1=I 
END  IF 

IF  (DGTEST(I) . GT.BPRIN2. AND. NSTOP . EQ. 1. ) THEN 
BPRIN2=DGTEST(I) 
IBIG2=I 
END  IF 

QTEST= ( 1 . - ABS ( REAL (Z ( I , ISMALL ) ) ) ) 
IF (QTEST .GE . BPRIN3 )THEN 
BPRIN3=QTEST 
IBIG3=I 
END  IF 
130      CONTINUE 
C 
C 

Qkkkkkkkkkkkkkkkk    OUTPUT  ****************** 
C 

IF (NSTOP. EQ.O.) THEN 

WRITE(8,900)IBIG1,BPRIN1,IBIG3,BPRIN3 
END  IF 
IF (NSTOP. EQ.1.)THEN 

WRITE (8 , 910 ) IBIG1 , BPRIN1 , IBIG2 , BPRIN2 , IBIG3 , BPRIN3 
END  IF 
C 

Qkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkk 

C*****  THE  FOLLOWING  THREE  LINES  OF  CODE  ARE  THE  ONES  WHERE  ***** 
Qkkkkk  THE  INTEGER  VARIABLE  LR  IS  ASSIGNED  THE  ORDER  OF  THE  ***** 
C*****  LEAST  OBSERVABLE  STATE.  IF  THE  USER  WANTS  TO  CHANGE  ***** 
C*****  OTHER  ENTRY  OF  C  RATHER  THAN  THE  DEFAULT  ONE  THESE  3  ***** 
Ckkkkk   LINES  SHOULD  BE  COMMENTED  AND  THE  4TH  ONE  ACTIVATED   ***** 

Qkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkk 
C 

IF(K.EQ.N)THEN 

LR=IBIG3 
END  IF 
C         LR=ORDER  OF  THE  STATE  TO  BE  CHANGED 


C 


TEST(IBIG1)=-1. 
DGTEST(IBIG2)=-1. 

Z(IBIG3,ISMALL)=CMPLX(-1. ,0.) 
140   CONTINUE 
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CALL  COMDM(QT,M,N,CNQ) 

WRITE(8,610TCNQ 

WRITE(6,705) 

705   FORMAT ('  '  '  FINISHING  THE  DEGREES  OF  OBSERVABILITY') 
IF(NT.EQ.l)GO  TO  999 
C 

c 

c*****  COMPUTE  THE  LEFT  EIGENVECTORS  OF  A  ****** 
C 

CALL  EIGRF(A,N,10,1,W,Z,10,WK1,IER) 
C 
C 

C*****  TRANSFORM  W  IN  P  WHERE  P  ARE  THE  REAL  PARTS  OF  W  ***** 
C 

DO  1000  1=1, N 

DO  990  J=1,N 

P(I,J)=REAL(Z(I,J)) 
990       CONTINUE 
1000  CONTINUE 
C 
C 

C*****  FIND  THE  MINIMUM  CONDITION  NUMBER  FOR  Q  TRANSPOSED  ***** 
C********  AND  CORRESPONDING  NEW  ENTRY  IN  THE  C  MATRIX  ********* 
C 

WRITE(6,710) 
710   FORMAT ('  ','  MINIMIZING  THE  CONDITION  NUMBER  OF  QT ' ) 

Y1=CNQ 

CALL  MULT2(C,P,N,G,AL,LS,LG) 

LSOLD=LS 

LGOLD=LG 

CALL  CONDN(G,N,N,CNG) 

X1=CNG 
C 

c 

C******  COMPUTE  THE  CNQ  FOR  A  CNG  THAT  EQUALS  ONE  ***** 
C 

PK=1. 

CALL  CNEW(C,P/LR,LS,LG,PK,N,CN) 

C(1,LR)=CN 

CALL  MULT1(C,A,N,QT) 

CALL  CONDN(QT,N,N,CNQ) 

Y2=CNQ 

X2=PK 
C 

280   CONTINUE 
C 

c 

C*****  FIND  THE  ORDER  OF  THE  LARGEST  AND  SMALLEST  ENTRY  IN  G  ***** 
C 

CALL  MULT2(C,P,N,G,AL,LS,LG) 
C 

c 

C**********  BISSECT  THE  INTERVAL  XI, X2   ********* 
C 

X3=(Xl+X2)/2. 

CNG=X3 
C 
C 

c*****  FIND  NEW  ORDINATE  (CNQ)  FOR  THE  PREVIOUS  ABCISSA  (CNG)  ***** 
C 

PK=1./CNG 

CALL  CNEW ( C , P , LR , LSOLD , LGOLD , PK , N , CM ) 

C(1,LR)=CN 

CALL  MULT1(C,A,N,QT) 

CALL  CONDN(QT,N,N,CNQ) 

Y3=CNQ 
C 
C 
C*****  COMPUTE  THE  SLOPE  NEAR  X3  ***** 
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XDEL=X3+DEL 
CNG=XDEL 
PK=1./CNG 

CALL  CNEW ( C , P , LR , LSOLD , LGOLD , PK , N , CN ) 
C(1,LR)=CN 
CALL  MULT1(C,A,N,QT) 
CALL  CONDN(QT,N,N,CNQ) 
YDEL=CNQ 

SLOPE=(Y3-YDEL)/(X3-XDEL) 
C 
C 

c*****  RENAME  THE  NEW  POINTS  SUCH  THAT  IN  THE  INTERVAL  ***** 
c*****  TOWARDS  THE  DIRECTION  OF  DESCENDING  SLOPE  ***** 
c*****  THE  point  WITH  THE  LARGEST  ORDINATE  IS  X1,Y1  ***** 
c*****  AND  THE  ONE  WITH  THE  SMALLEST  ORDINATE  IS  X2,Y2  ***** 
C 

IF(SLOPE.GE.O. )THEN 
IF(Y2.LT.Y1)THEN 

IF(Y3.LT.Y2)THEN 
Y1=Y2 
X1=X2 
Y2=Y3 
X2=X3 
ELSE 

Y1=Y3 
X1=X3 
END  IF 
ELSE 

IF(Y3.LT.Y1)THEN 
Y2=Y3 
X2=X3 
ELSE 

Y2=Y1 
X2=X1 
Y1=Y3 
X1=X3 
END  IF 
END  IF 
ELSE 

IF(X2.LT.X1)THEN 

IF(Y3.LT.Y1)THEN 
Y2=Y3 
X2=X3 
ELSE 

Y2=Y1 
X2=X1 
Y1=Y3 
X1=X3 
END  IF 
ELSE 

IF(Y3.LT.Y2)THEN 
Y1=Y2 
X1=X2 
Y2=Y3 
X2=X3 
ELSE 

Y1=Y3 
X1=X3 
END  IF 
END  IF 
END  IF 

IF(0.5*ABS(X2-X1) .GT.l.E-5)GO  TO  280 
WRITE (6, 720) 
720   FORMAT ('  ','  OPTIMIZATION  COMPLETED ■ ) 

WRITE(6,725) 
725   FORMAT ('  ','  COMPUTING  THE  NEW  DEGREES  OF  OBSERVABILITY1) 
C 
C 
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Qv.kkkkkkkkkkkkkk    f,J£^/  OUTUP  ******************* 
C 

WRITE(3,500) 

WRITE(8,350)(C(1,J) ,J=1,N) 
C 

WRITE (8, 600) 
L13=13 
NT=NT+1 

IF(NT.EQ.l)  GOTO  26 
C 
C 

q   kkkkk    FORMATS  ***** 
C 
350   FORMAT ('  '  , F13 . 7 , 7 (2X, F13  .  7 
370   FORMAT ('  ' , F13 . 7 , 7 (4X , F13 .7 

THE  ENTRIES  OF  MATRIX  C  ARE  : ' ) 

THE  NEW  ENTRIES  OF  MATRIX  C  ARE  : ' ) 

THE  ENTRIES  OF  MATRIX  A  ARE  : ' ) 

Q  TRANSPOSE  IS  : ' ) 

Q  TRANSPOSE  NORMALIZED  IS  : ' ) 

THE  NEW  DEGREES  OF  OBSERVABILITY  ARE  : ' ) 

THE  CONDITION  NUMBER  OF  Q  TRANSPOSED  IS  EQUAL  TO 


450  FORMAT ('  ',//, 

500  FORMAT ('  ',//, 

550  FORMAT ('  ',//, 

560  FORMAT ('  ',//, 

570  FORMAT {'  ',//, 

600  FORMAT ('  ',//, 

610  FORMAT ('  ',//, 

&'  ,F13.7) 

800  FORMAT ( '  ' ,//, 

850  FORMAT!1  ABLIN 


STATE  X', 12, IX,' IS  NOT  OBSERVABLE ',// ) 
,'S  CRITERION  CAN  NOT  BE  EVALUATED,  ',/,'  BECAUSE 
&Q  TRANSPOSE  IS  SINGULAR1,/) 
880   FORMAT ('  ',/,'  RELATIVE  DEGREE  OF  OBS .  FOR  GRAPHICAL  ABLIN1,1  AND 

&  GENERALIZED  METHODS') 
900   FORMAT ('  ' , 2X, ■ X ' , 12 , 2X, F9 . 6 , 2X, 12 ( ' . ' ) , 2X, ' X' , 12 , F9 .6) 
910   FORMAT ( '  ' , 2X, ' X ' , 12 , 2X, F9 . 6 , 2X, ' X ' , 12 , 2X, F9 . 6 , 2X, ■ X ■ , 12 , 2X, F9 . 6) 
999   WRITE(6,730) 

730   FORMAT ('  ','  PROGRAM  COMPLETED ' ) 
STOP 
END 
c 
c 

Qkkkkkkkkkkkkkkkkkk     SUBROUTINES  ******************** 
C 

Qkkkkkk   SUBROUTINE  TO  COMPUTE  THE  COLUMN  NORMS  ***** 
C 

Qkkkkkkkkkkkkkkkkkk    PURPOSE  ************************ 

C*   THE  PURPOSE  OF  THIS  SUBROUTINE  IS  TO  COMPUTE    * 
C*THE  MAXIMUM  SUM  OF  THE  COLUMNS  OF  A  GIVEN  MATRIX  * 

Qkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkk 
C 

SUBROUTINE  NORM (DA ,NL ,NC, DBIG) 
REAL*4  DA(10,10) 
DBIG=-1. 
DO  3  J=1,NC 
TOC=0. 
DO  1  1=1, NL 

TOC=TOC+ABS(DA(I,J)) 
1  CONTINUE 

IF(TOC.GT.DBIG)THEN 

DBIG=TOC 
END  IF 
3         CONTINUE 
RETURN 
END 
c 
c 

c*****  SUBROUTINE  MULT1  (COMPUTES  Q  TRANSPOSED)  **** 
C 

Qkkkkkkkkkkkkkkkkkk    PURPOSE  ************************ 

C*   THE  PURPOSE  OF  THIS  SUBROUTINE  IS  TO  COMPUTE    * 
C*THE  PRODUCT  OF  THE  C  AND  A  MATRICES  TO  PRODUCE  QT* 

Qkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkk 

C 

SUBROUTINE  MULT1 (C , A,N ,QT) 
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REALM  C(l,10)  ,A(10,10),QT(10,10) 
DO  30  J=1,N 

QT(1,J)=C(1,J) 
30    CONTINUE 
C 

DO  60  1  =  2, N 

DO  50  JL=1,N 

QT(I,JL)=0. 
DO  40  JC=1,N 

QT(I,JL)=QT(I,JL)+QT(I-1,JC)*A(JC,JL) 
40  CONTINUE 

50        CONTINUE 
60    CONTINUE 
RETURN 
END 
c 
c 

C****  SUBROUTINE  MULT2  (COMPUTES  THE  MATRIX  G)  AND  GIVES  THE  ORDER  **** 
c*********  0F  THE  LARGEST  AND  SMALLEST  ENTRY  IN  THE  G  MATRIX  ********** 
C 

£*******  ***********  PURPOSE  ************************ 

C*  THE  PURPOSE  OF  THIS  SUBROUTINE  IS  TO  COMPUTE  * 
C*THE  PRODUCT  OF  THE  C  AND  P  MATRICES  TO  PRODUCE  G  * 
C*  IT  ALSO  PROVIDES  THE  ORDER  OF  THE  SMALLEST  AND  * 
C*  LARGEST  ENTRIES  IN  THE  G  MATRIX  AS  WELL  AS  AL  * 
C*  WHICH  IS  A  ROW  VECTOR  WHOSE  ELEMENTS  ARE  THE  * 
C*  ENTRIES  OF  G  * 

£*************************************************** 
C 

SUBROUTINE  MULT2 (C,P ,N .G, AL , LS , LG) 
REALM  C(1,10),P(10,10),AL(1,10),G(10,10) 
DO  1090  JL=1,N 
AL(1, JL)=0. 
DO  1080  JC=1.N 

AL(1,JL)=AL(1,JL)+C(1,JC)*P(JC,JL) 
1080     CONTINUE 
1090  CONTINUE 

GRL=-l.E+4 

PQL=l.E+4 

DO  2020  J=1,N 

IF(AL(1,J) .GT.GRL)THEN 
GRL=AL(1,J) 
LG=J 
END  IF 

IF(AL(1,J) .LT. POL) THEN 
PQL=AL(1,J) 
LS=J 
END  IF 
2020  CONTINUE 

DO  1094  1=1, N 

DO  1092  J=1,N 

IF(I.EQ.J)THEN 

G(I,J)=AL(1,J) 
ELSE 

G(I,J)=0. 
END  IF 
1092     CONTINUE 
1094  CONTINUE 
RETURN 
END 
c 
c 

c****  SUBROUTINE  CONDN  (COMPUTES  THE  CONDITION  NUMBER  OF  AN  ARRAY)  **** 
C 

SUBROUTINE  CONDN  (DAR,NL ,NC , DCN) 
REALM  DAR(NL,NC)  ,DARI(10,10)  ,WK2(130) 
CALL  NORM(DAR,NL,NC,BG) 
CALL  LINV2F ( DAR ,NL , 10 , DARI , 2 , WK2 , IER) 
CALL  NORM ( DARI, NL,NC,BGI) 
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DCN=EG*BGI 
RETURN 
END 
C 

c 

c*****  SUBROUTINE  TO  COMPUTE  THE  NEW  ENTRY  IN  THE  C  MATRIX  ***** 
C 

SUBROUTINE  CNEW ( C , P , LR , LS , LG , PK , N , CN) 

REAL*4  C(1,10),P(10,10) 

ADD=0. 

DO  2040  1  =  1, N 

IF(I.NE.LR)THEN 

ADD=ADD+ ( PK*P ( I , LG ) -P ( I , LS ) ) *C ( 1 , I ) 
END  IF 
2040  CONTINUE 

CN=(-ADD)/(PK*P(LR,LG)-P(LR,LS)) 

RETURN 

END 
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